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FOREWORD 


The scope of this study was to plan, develop, 
and implement methods for analyzing and improving 
the performance of digital data smoothing filters 
used in the Radar Data Reduction System at the Wallops 
Flight Facility (WFF). During the study the following 
primary objectives were accomplished: (1) The accu- 

racy of the current WFF data smoothing technique was 
analyzed for a variety of radars and payloads, using 
tracking data provided by WFF for this purpose; 

(2) alternative data noise reduction techniques were 
assessed and recommendations were made for improving 
radar data processing at WFF; (3) a data-adaptive 
algorithm, based on Kalman filtering and smoothing 
techniques, was developed for estimating payload 
trajectories above the atmosphere from noisy time- 
varying radar data; (4) the new trajectory estima- 
tion algorithm was tested and verified using radar 
tracking data from Peru provided by WFF. 

Significant contributions to this study were 
made by the following individuals: A.R. Leschack 

provided the algorithm for computing Keplerian trajec- 
tories in Section 5.3 and provided an independent 
numerical check of the nominal- trajectory algorithm 
in Section 5.4; J.D. Goldstein provided the discrete- 
time linearized equations of motion for small per- 
turbations about the Keplerian trajectory used in 
Section 5.4; and A.E. Rhenals developed software for 
analyzing the WFF data tapes, performed a significant 
part of the radar data error analysis under Tasks 1 
and 2, and contributed theoretical analyses of the 
effects of smoothing filters' on nominal payload 
trajectory signals and data noise covariances. 
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INTRODUCTION 



1.1 OBJECTIVES OF THE STUDY 


The scope of this study was to plan, develop, and 
implement methods for analyzing and improving the performance 
of digital data smoothing filters. These filters are used in 
the Radar Data Reduction System at the Goddard Space Flight 
Center/Wallops Flight Facility (GSFC/WFF) to reduce noise levels 
in radar tracking data. The study had four primary objectives: 


• To develop stochastic models for radar 
tracking data provided by WFF 

• To determine the propagation of the radar 
data errors through the WFF noise-reduction 
filters into positional data products 

• To assess alternative noise reduction 
techniques and to make recommendations 
for improving the current filtering 
techniques 

• To develop and verify an algorithm for 
smoothing radar tracking data to estimate 
trajectories above the atmosphere. 


To meet these objectives, the study was divided into 
four tasks: 


• Task 1 - Analysis of Radar Data 

• Task 2 - Analysis of Existing Filters 
and Error Propagation into Positional 
Data 
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• Task 3 - Assessment of Alternative Noise- 
Reduction Techniques 

• Task 4 - Development and Verification of 
Trajectory Estimation Algorithm. 

1.2 TECHNICAL APPROACH 

The analysis of WFF radar data and smoothing filters 
(Tasks 1 and 2) is based on stochastic modeling and covariance 
analysis of noise-like errors. A combination of autoregressive 
and state-space techniques is used. The autoregressive models 
provide a cross-check on the more flexible state-space models. 

To estimate the root-mean-square (rms) noise levels, state-space 
covariance analysis is used. 

The assessment of alternative smoothing techniques 
(Task 3) is based on the results of the first two tasks and 
the theory of Kalman optimal filtering and smoothing. 

A new algorithm for processing noisy radar data for 
trajectories above the atmosphere is developed and verified 
under Task 4. The algorithm is based on Kalman filtering and 
smoothing techniques. Radar tracking data from Peru (provided 
by WFF) are used to verify the effectiveness of the algorithm 
with real data, including data from two radars simultaneously 
tracking a single payload. 

1.3 ORGANIZATION OF REPORT 

Each task of the study is documented, in a separate 
chapter (Chapters 2-5). The report ends with Chapter 6, which 



provides a summary of the investigation, the primary conclu- 
sions, and recommendations for further study. Supporting 
mathematical definitions and analyses are presented in Ap- 
pendices A, B, C, D, and E. In particular, the coordinate 
systems used in this study are defined in Appendix D. 
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TASK 1 - ANALYSIS OF WFF RADAR DATA 


2 . 1 INTRODUCTION 

2.1.1 Objective 

The main objective of Task 1 is to develop error models 
for WFF radar tracking data . The data consist of azimuth, 
elevation, and range measurements for a variety of payloads 
and radars A3 provided by the WFF PASS-1 data processing program 
(Ref. 1). The noise-like zero-mean errors in the radar data 
are modeled as stochastic processes. In contrast, the system- 
atic errors in the radar data and the payload trajectories are 
modeled as polynomials. 

2.1.2 Radar Data 

The radar data sets analyzed in this study are listed 
in Table 2.1-1. The first three radars listed in the table 
(Radars Nos. 3, 5, and 6) correspond to data analyzed under 
Task 1. The remainder of the data sets (Radars Nos. 8 and 41) 
are analyzed under Task 4. 

2.1.3 Approach 

The technical approach for analyzing the radar data 
consists of three steps: 

• Model the nominal trajectory for each 

time series using least-squares orthogonal 
polynomial functions of time 
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TABLE 2.1-1 

RADAR DATA SETS ANALYZED 


RADAR 

RADAR 
LOCATION 
N . LATITUDE 
E. LONGITUDE 
HEIGHT 

TRAJECTORY 

NAME , MODEL NO . , DATE 

NO. 3 
(WFF) 

37.841309 deg 
-75.485102 deg 
14.08 m 

ZUN1, El-0425... 0427, 12/1/82 

S. LOKI OPTICAL, (MODEL & DATE UNKNOWN) 

NO. 5 
(WFF) 

37.860229 deg 
-75.509309 deg 
16.66 id 

S. LOKI SPHERE, Tl-0503, 4/25/83 
S. LOKI OPTICAL, (MODEL & DATE UNKNOWN) 

NO. 6 
(WFF) 

37.841585 deg 
-75.484692 deg 
9.43 m . 

S. LOKI OPTICAL, Tl-6615, 12/11/81 

NO. 8 
(PERU) 

-12 . 4993 deg 
-76.7965 deg 
74.26 m 

NIKE-ORION 31.027, TU2-0367, 3/9/83 
TERRI ER-MALEMUT£ 29.019, TU-0364, 3/15/83 

NO. 41 
(PERU) 

-12.4990 deg 
-76.7954 deg 
71.02 m 

NIKE-ORION 31.027, TU2-0367 , 3/9/83 


• Select the appropriate degree for each 
polynomial by using a spectral analysis 
of the residual tracking data (residual = 
data - polynomial) 

• Develop stochastic state-space models 
for the noise-like residual tracking 
data (azimuth, elevation, and range). 


Detailed information about the technical approach is presented 
in Section 2.2. 


2.1.4 Interpretation 

The polynomials fitted to the data are nominal estimates 
of the payload trajectory plus any systematic measurement errors. 
The residual data (actual measurements minus a polynomial) are 
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a combination of measurement errors and zero-mean noise-like 
signals caused by payload dynamics. The statistics of the 
residual data are used to model the noise in the radar measure- 
ments. 


2.2 TECHNICAL APPROACH 

2.2.1 Data Signal Components 

The radar tracking data may be represented as the sum 
of three signal components: 

• 

Tracking _ Payload + Systematic Noise-Like 

Data Motion Error Error 

( 2 . 2 - 1 ) 

The goal of the analysis is to estimate the statistics of the 
noise-like errors in real tracking data . This requires that 
the noise-like error signals be distinguished from the payload 
motion and systematic error signals. The following criteria 
are used to distinguish between the data signals: 


Payload Motion 

( • 

Polynomial- like 

and 

1 


Systematic Error 

[• 

High Coherence between Azimuth, 
Elevation, and Range Data 


| • 

Zero-Mean Random Noise 

Noise-Like Error 

!. 

Low Coherence between Azimuth, 
Elevation, and Range Data 


These criteria are the logical basis for the data analysis 
described in the next section. 
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2.2.2 Data Analysis 


The radar tracking data are analyzed in three steps: 


• Step 1 - Fit orthogonal polynomials to 
the raw tracking data . The purpose of 
this is to estimate the signal components 
caused by payload motion and systematic 
error. Orthogonal polynomials are used 
to avoid numerical problems that can 
otherwise make it difficult to fit the 
polynomials accurately. 

• Step 2 - Develop stochastic models for the 
residual tracking data . The residual 
data ( raw data minus polynomial) are 
expected to consist mostly of noise-like 
radar measurement errors when the degree 
of the polynomial is appropriate. The - 
stochastic models for the residuals are 
developed using autoregressive and state- 
space modeling techniques. 

• Step 3 - Select appropriate polynomial 
degrees based on power spectra and spectral 
coherences . The stochastic models developed 
in Step 2 are used to estimate the spectra 
and coherences. The appropriate polynomial 
degree is selected so that the azimuth, 
elevation, and range residuals have small 
coherences and nearly flat power spectra 

at low frequencies. 


A block diagram of the data processing is shown in 
Fig. 2.2-1. Each block in this figure is explained in the 
following discussion. 

Inputs - In Fig. 2.2-1 the inputs are "raw data." 

This means that the input data are in the form provided by the 
PASS-1 data processing program used at WFF. After plotting 
the individual time series for azimuth, elevation, and range 
data, obviously incorrect data are replaced by reasonable values 
using simple linear interpolation. If a low-frequency analysis 
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Figure 2.2-1 Block Diagram of Data Analysis 

is being performed, the optional resampling is performed at a 
slower sampling rate than the raw data. No anti-aliasing filter 
is used in this resampling operation. 

Fitting Polynomials - Orthogonal polynomials are fitted 
to subsets of the resampled data using least squares. The 
outputs of this procedure are: (1) the coefficients of the 

polynomial and (2) the residual data. The residuals are defined 
as the resampled data minus the polynomial. Typical subsets 
contain 500 measurements, and polynomial degrees usually range 
from 4 to 10 for the data sets analyzed in this study. The 
mathematical details of the polynomial fitting are presented 
in Appendix A. 

Modeling Residual Data - The residual data are analyzed 
using autoregressive and state-space stochastic modeling tech- 
niques. (Mathematical discussions of these techniques are 
presented in Appendices B and C.) The stochastic models are 
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shaping filters, i.e., difference equations driven by white 
noise. The steady-state solutions of these difference equa- 
tions define random processes that are models of the residual 
radar data. Power spectra and spectral coherences for the 
models are computed using formulas presented in Appendix B 
and Section 2.3.3. 

Polynomial Degree Selection - The power spectra and 
spectral coherences are examined to determine if the polynomial 
degree is too small or too large. If there is significant 
coherence (>25%) at low frequencies, then the degree is judged 
to be too small. And if the power spectrum is scooped-out at 
low frequencies, the degree is judged to be too large. T-hese 
criteria are consistent with the goal of distinguishing the 
signals caused by payload motion and systematic error from the 
signal caused by random measurement noise. As indicated in 
Fig. 2.2-1, if the degree needs to be changed, then the analysis 
process is repeated starting with a different polynomial degree. 
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2.2.3 Example Results 


As an example of the results obtained, Figs. 2.2-2 to 
2.2-4 depict the raw tracking data and residual data for Radar 
No. 3 tracking a Super Loki Optical payload. The ejection of 
the payload occurs at approximately 120 s. 

Azimuth Data - Figure 2.2-2 shows that azimuth was 
nearly constant at 134 degrees (deg). A degree-10 polynomial 
is appropriate for modeling the payload motion and systematic 
error over this segment of data. The residuals about this 
polynomial look like homogeneous random noise. There is no 
visible anomally caused by the payload ejection. 


Elevation Data - Figure 2.2-3 depicts the Loki-Optical 
elevation data. In this case a degree-8 polynomial was found 
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AZIMUTH DATA 



TIME Isecl 


Figure 2.2-2 


ELEVATION DATA 



TIME (seel 


Figure 2.2-3 


ELEVATION (deg) 


RESIDUAL AZIMUTH DATA 





Loki Optical Azimuth Data 


1 




RESIDUAL ELEVATION DATA 


. o 



,oki Optical Elevation Data 
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RANGE (kft) 


to be appropriate. The residuals appear to be more random 
than the azimuth residuals shown in Fig. 2.2-2, which was con- 
firmed by spectrum analysis. (The more random appearance 
corresponds to more power at higher frequencies.) 

Range Data - The Loki-Optical range data are shown in 
Fig. 2.2-4. In this case a degree-8 polynomial was appropriate. 
Although the payload ejection is not visible in the raw range 
data, a strong localized inhomogeneity caused by the ejection 
is seen in the residual data. This is an example where the 
residual data contain a combination of both measurement noise 
and residual payload motion. 


RANGE DATA RESIDUAL RANGE DATA 


R ORfiS*. 



TIME (sec) TIME (sec) 


s 


V 


Figure 2.2-4 Loki Optical Range Data 
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Loki-Sphere Range Data - Figure 2.2-5 shows data from 
Radar No. 5 tracking a Super Loki Sphere payload. Two sets of 
analysis results are included for comparison. The pair on the 
left in Fig. 2.2-5 shows the residual range data for a degree-6 


polynomial and the power spectrum estimated from the residual 
data using the autoregressive modeling technique described in 


Appendix B. The polynomial degree in this case is not too 
large because it yields a nearly flat power spectrum at low 
frequencies. In contrast, the right pair of plots show the 
results of using a polynomial of degree 15, which is too large . 
The degree-15 polynomial partially fits the low-frequency ran- 
domness of the radar data, as indicated by the dip in the power 




RESIDUAL RANGE DATA RESIDUAL RANGE DATA 



37 93 37 93 

TIME (sec) TIME (sec) 


RESIDUAL POWER SPECTRUM RESIDUAL POWER SPECTRUM 



FREQUENCY (Hz) FREQUENCY (Hz) 


Figure 2.2-5 Example of Selecting the Trajectory 

Polynomial, Loki Sphere Tracking Data 
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spectrum. The result of using the degree-15 polynomial would 
be to slightly underestimate the low-frequency noise in the 
radar data. Therefore, the degree-6 polynomial is preferred. 


2.3 ANALYSIS RESULTS 

2.3.1 Introduction 

Section 2.3 presents the results of applying the data 
analysis described in Section 2.2 to a variety of data sets 
provided by WFF. The principal results of this analysis are 
stochastic state-space models for residual radar tracking data . 
These models are used to estimate the power spectra, spectral 
coherences, and covariance matrices of noise-like errors in the 
tracking data. The covariance matrices are used in Chapter 3 
to estimate the rms random errors in positional data products 
expressed in latitude, longitude, and height. 

2.3.2 State-Space Models 

In this secton the concept of a state-space stochastic 
model is introduced. Mathematical details about the state-space 
modeling technique used in this study are provided in Appendix C. 


State-Space Equations - The residual tracking data 
consist of three time series, one for each of the azimuth, 
elevation, and range measurements. The three channels of data 
are combined to form a 3-vector: 


*k = 


azimuth 

elevation 

range 


at time k = 1,2,3... 


(2.3-1) 
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A state-space model for the vector random process y k is repre- 
sented by two equations: (1) a vector difference equation for 

the states; and (2) an algebraic equation relating the states 
to the observed process y k . The algebraic equation is: 

*k = ^k + 2 k 

x k = n*l matrix of n state varibles 
H = 3*n output matrix 
£ k = 3x1 matrix of 3 white-noise processes 

Cov(v^) = R = 3 X 3 covariance matrix (2.3-2) 

The white noise is called the innovations vector because it 
represents that part of the residual radar data y k which is 
uncorrelated with the past radar data; the term Hx^ in Eq. 2.3-2 
represents that part of y k which is correlated with the past . 
The state vector x k contains all necessary information about 
the past data (z^-i » Zfc-2 >•••)• 

The state vector satisfies the following difference 

equation: 

<frx k + Gv k . 

nxn transition matrix 

nx3 noise-gain matrix (2.3-3) 

Equations 2.3-2 and 2.3-3 represent a state-space model in the 
innovations form. There are other forms of state-space models, 
but this is the form that is appropriate for stochastic modeling 
based on empirical data. 


-k+1 = 
G 
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Modeling - The method of stochastic state-space mod- 
eling used in this study is based on a canonical-variates (CV) 
analysis of the data. (Mathematical details are described in 
Appendix C.) The canonical variates are used as estimates of 
the state variables: 




1st canonical variate 
2nd canonical variate 


nth canonical variate 


most important state 
2nd most important 


least important state 
(2.3-4) 


By defining the state variables with the most important state 
listed first and the least important state listed last, it is 
straight-forward to compute a family of models by adding one 
state (canonical variate) at a time. In this way, a CV anal- 
ysis of the radar data that yields n canonical variates can be 
used to compute the 4> , G, H, and R matrices for n+1 different 
models. Each model contains a different number of state vari- 
ables, ranging from zero states (a pure white-noise model) to 
an n-state model of maximum complexity. The Akaike information 
criterion (AIC) (Refs. 2-4) is then used to select from this 
family the one model that is best justified by the finite 
amount of radar data that was used for the analysis. This 
procedure avoids the problems of under modeling with too few 
states or over modeling with too many states. 

2.3.3 Applications of State-Space Models 


Power spectra - In this study, the state-space models 
are used to estimate power spectra, spectral coherences, and 
covariances of residual radar data. The power spectral density 
matrix of the residual radar data is a 3x3 matrix: 
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S AA S AE S AR 
S yy. = S EA S EE S ER 
_ S RA S RE S RR 


A ~ AZIMUTH 


E ~ ELEVATION 


R ~ RANGE 


(2.3-5) 


The elements along the main diagonal (S^, Sgj., and S RR ) are the 
auto spectra for the azimuth, elevation, and range data. The 
off-diagonal elements are cross spectra between the indicated 
pairs of measurements. The spectral density matrix is computed 
as follows using the state-space parameter matrices , G, H, 
and R: 


S yy (F) = H[ Ie i2nF - <J> ] 


1 

G+IJ R 


I+G T [Ie“ l27tF - <|> T ] 


1 

H 1 


(2.3-6) 


F = Normalized Frequency [cycle/sampling interval] 
I = Identity Matrix (either n*n or 3x3) 


The elements of matrix S yy are expressed in the same units as 
the corresponding elements of the covariance matrix of The 

normalized frequency F ranges from -1/2 to 1/2 [cycle/sampling 
interval] and is related to the sampling frequency f sam p [Hz] 
and the regular frequency variable f [Hz] as follows: 


F = f/f 


samp 


(2.3-7) 


The spectrum S yy (F) has the units of [variance]. It can be 
scaled to have the units of [variance/Hz], which are convention- 
ally used for continuous- time signals, and then expressed as 
S yy (f), a function of the frequency variable f [Hz] as follows: 


S (f) = 

yy 


A f/f. 


samp 


(2.3-8) 
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The covariance matrix P of y, is equal to the integral of 

yy K 

the power spectrum (integration of the spectral density matrix 
element-by-element): 


■ / 


S yy (D dF = 


/samp / 2 
Syy(f) ^ 

' f samp / 2 


(2.3-9) 


Spectral Coherence - The cross spectra (off-diagonal 
elements of S ) are complex valued because they represent the 
phase of the cross correlations between pairs of measurement 
processes. To suppress the phase information and focus attention 
on the magnitudes of the correlations as a function of frequency, 
the squared coherence (spectral coherence) is computed for 
each pair of measurement processes. For example, the squared 
coherence between the azimuth and range measurements, 
is defined by the formula: 


C AR< F) - 


' s ar< f >I 


(2.3-10) 


The coherence between any other pair of measurements is defined 
in the same way. 

The spectral coherence ranges in value from 07o to 
100%. It measures how much of the variance of one random proc- 
ess can be explained as a linear transformation of another 
random process. Put differently, the coherence measures how 
much power the two processes have in common with each other on 
the average. 


The coherence is a function of frequency and provides 
more information about the crosscorrelation structure than 
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simple correlation coefficients. (The squared correlation 
coefficient is equal to the area under the squared coherence 
function.) The spectral information contained in the coherence 
function is used in this study to determine when residual radar 
data from two different measurements (e.g., azimuth and range) 
contain correlated low-frequency signals (polynomial-like tra- 
jectory signals or systematic error signals). A significant 
occurrance of these signals produces coherence larger than 25% 
and indicates that a higher-degree polynomial is appropriate 
for computing the residual data. 

Covariance Matrices - The state-space models are also 
used to compute covariance matrices for the residual radar 
data. The method for computing covariances described in this 
section is much more convenient than evaluating one of the 
integrals in Eq. 2.3-9. These covariances are used in Task 2 
to estimate the rms errors in positional data derived from the 
radar measurements. 

The covariance matrix of the residual radar measure- 
ments is denoted P . It is computed by solving the following 
matrix equations, which contain the state-space parameter 
matrices <p , G, H, and R. 


P 

yy 


cov <x k ) 


Covariance Matrix of Azimuth, 
Elevation, and Range Data 


[3x3] 


P xx = cov(x^) = Covariance Matrix of State Vector [nxn] 


P yy = H P xx H? + R (2.3-11) 

P xx = ^ P xx * T + GRgT (2.3-12) 


A practical way of solving equation 2.3-12 is to set 
P xx initially equal to the nxn identity matrix and use this 
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matrix to evaluate the right side of Eq. 2.3-12. Then the new 
value of P is used to evaluate the right side again, etc. 
until the value of remains unchanged to within the desired 
computational accuracy. This algorithm converges if <|> repre- 
sents a stable state-space model (i.e., all eigenvalues of <j> 

have moduli less than unity). Once P is computed, P is 

xx yy 

computed using Eq. 2.3-11. 

2.3.4 Error Models 

In this section state-space models are presented for 
noise-like radar measurement errors. These models are based 
on analyses -of residual radar data from WFF Radars Nos. 3, 5, 
and 6, which were tracking Zuni , Loki Optical, and Loki Sphere 
payloads. Root-mean-square (rms) values, power spectra, and 
spectral coherences are discussed. 

Estimated rms Error Levels - Table 2.3-1 summarizes 
the rms values of the residual tracking data from four data 
sets. These rms values are estimates of the rms noise levels 
in the tracking data. The rms values were computed from 
stochastic state-space models for residual tracking data and 
do not include systematic errors in the tracking data . The 
residual angle data from radars Nos. 3 and 5 have rms values 
ranging from 2.9 mdeg to 11.4 mdeg, while the rms range data 
are more tightly clustered from 5.5 ft to 6.8 ft. Radar No. 6 
has significantly higher estimated noise levels because it is 
a wide-band acquisition radar, which is not intended for accu- 
rate tracking. The Zuni error model is based on the analysis 
of three trajectories, while the other models are based on 
segments from single trajectories. More information about the 
estimated noise levels is presented in the following discussions 
of power spectra and spectral coherences. 
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TABLE 2.3-1 

ESTIMATED RMS NOISE LEVELS 



NOMINAL 

RMS 

DATA SET 

PAYLOAD 

RESIDUAL 


COORDINATES 

TRACKING DATA 




ZUNI 



(Radar #3) 



(3 Trajectories) 



Azimuth 

132 deg 

4.4 mdeg 

Elevation 

15 deg 

6 . 2 mdeg 

Range 

27 kft 

5.5 ft 

LOK1 SPHERE 



(Radar #5) 



(37s to 93s) 



Azimuth 

137 deg 

2 . 9 mdeg 

Elevation 

42 deg 

4 . 7 mdeg 

Range 

27 kft 

6.0 ft 

LOKI OPTICAL 



(Radar #3) 



(80s to 130s) 



Azimuth 

139 deg 

11.4 mdeg 

Elevation 

78 deg 

11.4 mdeg 

Range 

220 kft 

6.8 ft 

LOKI OPTICAL 



(Radar #6) 



(22-27 min) 



(33-38 min) 



Azimuth 

98 deg 

53 mdeg 

Elevation 

18 deg 

41 mdeg 

Range 

260 kft 

41 ft 


Radar No. 3 Power Spectra and Coherences - Figure 2.3-1 
depicts the estimated power spectra (PSDs or power ppectral 
densities) and squared coherences for the errors in the Zuni data 
from radar No. 3. These graphs were computed from the stochastic 
model derived from an analysis of the residual tracking data 
for three trajectories at a sampling rate of 10 samples per 
second. The coherence plots on the right side of Fig. 2.3-1 
indicate that the squared coherences are less than 25% at low 
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frequencies. At the same time, the spectrum plots on the left 
side are flat at low frequencies. This verifies that the low- 
frequency signals caused by payload motion and systematic 
measurement, errors have been appropriately modeled and removed 
from the residual data. 

According to Fig. 2.3-1 the range errors have the 
most nearly white (flat) spectrum, while the elevation errors 
are more nearly band limited. The frequency at which the 
elevation PSD decreases to half its low-frequency value is 
1 Hz. The corresponding half-power frequency for the azimuth 
errors is 2.5 Hz. 

Figure 2.3-2 depicts another set of estimated PSDs 
and coherences for noise-like errors in the data from radar 
No. 3. In this case, the radar was tracking a Super Loki 
rocket prior to ejection and its Optical payload after ejec- 
tion. These plots were computed using a 10-sample/second 
tracking data shown in Figs. 2.2-2 through 2.2-4. The local- 
ized inhomogeniety in the range data (Fig. 2.2-4) is caused by 
the payload ejection and produces the hump in the range PSD 
(Fig. 2.3-2) at 0.7 Hz. Aside from this anomaly, the range 
PSD in Fig. 2.3-2 for the Optical payload is similar to the 
range PSD in Fig. 2.3-1 for the Zuni trajectories. 

The azimuth and elevation PSDs in Fig. 2.3-2 for the 
Loki Optical tracking errors are significantly different from 
their Zuni counterparts in Fig. 2.3-1. The Loki Optical error 
spectra have 10-100 times as much power at low frequencies; 
while at high frequencies the two sets of PSDs are more nearly 
equal . 


Radar No. 5 Power Spectra and Coherences - Figure 2.3-3 
depicts the estimated PSDs and coherences for noise-like errors 
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in tracking data from radar No. 5, which was tracking a Super 
Loki Sphere payload. The state-space model used in this anal- 
ysis was developed from a segment of the tracking data extending 
from 37 s to 92.7 s with a sampling rate of 10 samples per 
second. The small coherences and flat PSDs at low frequencies 
indicate that polynomials of appropriate degree were used in 
computing the residual data. The flat portions of the azimuth 
and elevation PSDs at high frequencies indicate definite white- 
noise floors in the angular measurements. The magnitudes of 
these floors denote less high-frequency noise in Radar No. 5 
data as compared to corresponding noise levels in Radar No. 3 
data represented by Figs. 2.3-1 and 2.3-2. At low frequencies 
the Radar No. 5 data have noise power levels that are between 
those estimated for the Loki Optical and Zuni data from Radar 
No. 3. 


2.4 SUMMARY AND CONCLUSIONS 
2.4.1 Summary 

Under Task 1 of the study, the noise-like errors in 
radar tracking data were modeled using state-space techniques. 

The data sets included tracking radars Nos. 3 and 5, and acqui- 
sition radar No. 6. The trajectories analyzed were for Zuni, 
Super Loki Optical, and Super Loki Sphere payloads. 

Random measurement noise signals were separated from 
trajectory signals and systematic tracking errors in the data 
by subtracting least-squares orthogonal polynomials from segments 
of the data. The resulting residual data were used as estimates 
of the noise-like error signals in the tracking data. The 
appropriate degrees of the polynomials were determined from 
spectral analyses of the residual data using an autoregressive 
modeling technique. 
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2.4.2 Conclusions 


The analysis of data from tracking Radars Nos. 3 and 5 
and acquisition Radar No. 6, leads to the following conclusions 


• The estimated rms noise levels in tracking 
data from Radars Nos. 3 and 5 vary from 
2.9 millidegree (mdeg) to 11 mdeg in 
azimuth and elevation, and 5.5 ft to 
6.8 ft in range, depending on the radar 
and the payload (Zuni, Super Loki Optical, 
and Super Loki Sphere) 

• The estimated rms noise levels in Super 
Loki Optical data from the wide-band 
acquisition Radar No. 6 are 53 mdeg for 
azimuth, 41 mdeg for elevation, and 
41 ft for range. 
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3 . 


TASK 2 - ANALYSIS OF ERROR PROPAGATION 
INTO POSITIONAL DATA 


3 . 1 INTRODUCTION 

3.1.1 Objective 

The main objective of Task 2 is to determine the propa- 
gation of the radar data errors through the WFF noise-reduction 
filters into positional data products expressed in latitude, 
longitude, and height . The noise-reduction filters currently 
used at WFF are digital finite-impulse-response (FIR) smoothing 
filters (low-pass filters with symmetric impulse responses 
that produce zero phase shift). The parameters of the filters 
are manually selected to attenuate the high-frequency noise in 
the radar tracking data. The radar data are passed through 
these filters in the SMAD data processing program (Ref. 1), 
which produces smoothed tracking data as output. These smoothed 
data are then used as inputs to the MESUP and POSDAT programs 
(Ref. 1) that generate positional data products. The objective 
of this task is to estimate the rms noise-like errors in these 
data products, based on the data error models developed under 
Task 1. 

3.1.2 Approach 

The technical approach used in the analysis of error 
propagation consists of three steps, which are depicted in 
Fig. 3.1-1 and discussed in the following: 

• Step 1 - Develop state-space equations 

for the WFF smoothing filter being analyzed . 

The input to this step is the impulse 
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response of the WFF filter. The output 
is a set of state-space equations that 
are mathematically equivalent to the 
specified impulse response. 

• Step 2 - Compute the error covariance 
matrices of the smoothed radar tracking 
data . The inputs to this analysis are 
the equations from Step 1 and the stochas-. 
tic model for residual tracking data, 
which was developed under Task 1. The 
output is the covariance matrix for the 
noise-like errors in the smoothed residual 
tracking data. 

• Step 3 - Compute the error covariance 
matrices of positional data products 
expressed in latitude, longitude, and 
height . The inputs to this step are the 
covariances from Step 2, the radar posi- 
tion coordinates, and the nominal payload 
coordinates. The output is the error 
covariance matrix of the positional data. 



Figure 3.1-1 Block Diagram of Error Propagation Analysis 
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In addition to the analysis of noise-like errors, a 
systematic filtering error is also analyzed. This systematic 
error is caused by the slight distortions suffered by each 
payload trajectory signal as it is processed by the WFF filters 
This error is called filter- induced trajectory bias . It is 
analyzed in this study by computing the distortions suffered 
* by polynomials (which model nominal payload trajectory signals) 

when the polynomials are smoothed by the WFF filters. 

3.2 STATE-SPACE FILTER EQUATIONS 

The purpose of this section is to describe a state- 
space representation for WFF smoothing filters. This repre- 
sentation is convenient for computing the error covariance 
matrices of smoothed radar tracking data, given state-space 
models for the noise in the unsmoothed tracking data. > 

The impulse response h^ of a FIR filter is represented 
by the sequence of numbers 

h_^ > 2 • •••* ^-i’ h o> ^1* * * * * 1. * (3.2 — 1) 

In Eq. 3.2-1, L denotes the half-length of the impulse response 
The smoothed output y^ (for k = L+l , L+2 , ...) produced by an 
input sequence u^ (for k = 1, 2, ...) may be represented by 
the following convolution: 

L 

y k = E h j u u-j < 3 - 2 - 2 > 

j = -L 

An equivalent representation, and one that is more convenient 
than Eq. 3.2-2 for covariance calculations, is the following 

t 

state-space model with the n’xl state vector (the number of 
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states n' equals the number of samples in the support of the 
filter's impulse response sequence, i.e., n' = 2*L+1): 


Vu = H'x, 


-k+1 


= 0'x, + G'u, 


0* = 


0 1 0 0 

0 0 1 0 0 

0 0 0 1 0. . .0 


G' = 


0 

0-1 


0 


l lxl ] (3.2-3) 

fn'xl] (3.2-4) 


[n 'xn ' J 


(3.2-5) 


[n'xl] 


(3.2-6) 


H' = lh. L h _ L+ i ...hg...h L ] { l x n ' ] 


(3.2-7) 


Equation 3.2-3 states that the filter output is a linear 
function of the state vector x^; H' is a lxn' matrix, which 
contains the impulse response of the filter as defined by 
Eq. 3.2-7. Equation 3.2-4 is the state propagation equation 
in which 0' is the n'xn' transition matrix and G' is the n'xl 
input matrix. 

In Eqs. 3.2-2 and 3.2-4, the input sequence u^ repre- 
sents the noise-like error signal in radar tracking data (azimuth, 
elevation, or range measurements). Stochastic state-space models 
for u^ were developed under Task 1 as discussed in Chapter 2. 

Such noise models are represented here by the following equations 

IT 

involving the n"xl state vector x^: 
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u, = H"x, + v, , variance(v.) = R" [ lx 1 ) (3.2-8) 


*k+1 = ♦”£ 1 , + G " v i 


[n’’xl] 


(3.2-9) 


Equations 3.2-8 and 3.2-9 are a stochastic model for the radar 

tl 

noise u^; the scalar white noise v^ driving this model is called 

the innovations and is uncorrelated with x. for j < k. 

-J = 


For covariance calculations, the two state-space models 
for the WFF smoothing filter and the radar noise are combined 
to form one larger model of the following form, where n = n’+n" : 


*k+l 

= 0* k + G« k 

[nxl ] 

(3.2-10) 

*k 

= Hx k 

[lxl] 

(3.2-11) 


The matrices in Eqs. 3.2-10 and 3.2-11 are defined as follows: 



H = (H' 0] 


[nxl ] 

(nxn ] 
[lxn] 



(3.2-12) 

(3.2-13) 

(3.2-14) 

(3.2-15) 

(3.2-16) 
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3.3 ERROR COVARIANCE EQUATIONS 

The steady-state error variance of the smoothing filter 
output is denoted var(y). From a covariance analysis of 
the Eqs . 3.2-10 through 3.2-16 (Ref. 5), it can be shown that 
the variance of the noise in the filter output may be computed 
using the following equations: 

var(y) = H’P'H' 7 [lxl] ' (3.3-1) 

P' = 0'P'4>' T + <j>’SH" T G' T + G'H"S T 4»' T 

+ (G'H")P"(G’H") T + G'R"G ,T [n'xn’l (3.3-2) 

In Eq. 3.3-2, P' is the steady-state error covariance of the 

1 

filter state vector x^. To solve Eq. 3.3-2 for P', the matrices 
S and P" are first computed using the following equations: 

P" = <|)"P"<|>" T + G"R"G" T [n"x n " ] (3.3-3) 

S = 4» ’ S4> ,tT + G'H"P"4>" T + G'R"G" T [n'x n "] (3.3-4) 

The recommended procedure for solving these equations 
is to (1) solve Eq. 3.3-3 for P" , (2) solve Eq. 3.3-4 for S, 

(3) solve Eq. 3.3-2 for P’, and (4) use Eq. 3.3-1 to compute 
var(y). Equations 3.3-2 through 3.3-4 are equilibrium equa- 
tions describing statistical steady-state error covariances. 
They may be solved by using iteration, e.g., Eq . 3.3-3 may be 
solved by initially setting P" equal to the n"xn" identity 
matrix, evaluating the right side of the equation, using the 
new value for P" to re-evaluate the right side, and continuing 
until the elements of P" remain unchanged to within the desired 
numerical accuracy. This iterative method converges to the 
unique solution of the equation whenever the steady-state covar- 
iance exists (i.e., when the eigenvalues of P" have moduli 
less than unity). 
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3.4 


COORDINATE TRANSFORMATION EQUATIONS 


The last step in determining the rms errors in posi- 
tional data products is to transform the error variances of 
the smoothed tracking data (azimuth, elevation, and range) 
into the corresponding error covariances of the positional 
data (latitude, longitude, and height). The equations for 
this transformation are discussed in this section, while the 
mathematical details are presented in Appendix D. 


The coordinate transformation is performed in three 
steps as discussed in the following: 


• Step 1 - Transform from radar az-el-range 
coordinates to topographic north-east- 
down (NED) coordinates ! This is a non- 
linear transformation because radar 
coordinates are not Cartesian. 

• Step 2 - Transform from topographic NED 
coordinates to geocentric Cartesian 
coordiantes . This is a linear trans- 
formation from one Cartesian system to 
another . 

• Step 3 - Transform from geocentric 
Cartesian coordinates to geodetic lat-long- 
height coordinates . This is a nonlinear 
transformation . 


The propagation of the error covariance through this 
sequence of transformations is accurately approximated by 
linearizing the nonlinear transformations (in Steps 1 and 3) 
about the nominal payload coordinates. This approximation is 
accurate because the rms tracking errors are a small fraction 
of the nominal payload coordinates. 
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3.5 


ANALYSIS RESULTS 


This section presents rms error estimates for posi- 
tional data products expressed in latitude, longitude, and 
height. The error estimates were computed using the analysis 
techniques described in Sections 3.2 - 3.4 and the stochastic 
error models developed under Task 1. The accuracy estimates 
are for noise-like errors and apply to WFF radars Nos. 3 and 5 
tracking Zuni , Super Loki Sphere, and Super Loki Optical 
payloads. 

3.5.1 Rms Errors of Smoothed Tracking Data 

Rms error estimates for smoothed radar data (azimuth, 
elevation, and range measurements) are presented in Table 3.5-1. 
The first column indicates the payload and data sets used for 
analysis. (Data from wide-band radar No. 6 are excluded from 
this comparison because this radar is not intended for precise 
tracking applications.) The third and fourth columns give the 
estimated rms noise levels in the data before and after proc- 
essing with the WFF smoothing filter. (For these data sets, 
which were analyzed at a sampling rate of 10 samples/second, 
the appropriate WFF smoothing filter is designated by the code 
F00.040.10.) The last column in Table 3.5-1 indicates the 
percent reduction of rms noise in the smoothed data as compared 
with the rms noise before smoothing. 

The data in Table 3.5-1 lead to the following con- 
clusions. The percent reduction in estimated rms noise due to 
smoothing falls in the range of 71% to 14%. The estimated rms 
angular errors for Radar No. 3 vary significantly with payload: 
the Zuni trajectories have rms accuracies of 1.5 mdeg in azimuth 
and 2.7 mdeg in elevation; in contrast, the Loki Optical payload 
was tracked with accuracies of 9.5 mdeg in azimuth and 8.2 mdeg 
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TABLE 3.5-1 

STANDARD DEVIATIONS OF SMOOTHED RESIDUAL 
TRACKING DATA 


DATA SET 

NOMINAL 

PAYLOAD 

COORDINATES 

RMS 

BEFORE 

SMOOTHING 

RMS 

AFTER 

SMOOTHING 

PERCENT 
REDUCTION 
OF RMS 

ZUNI 

(RADAR #3) 

(3 TRAJECTORIES) 
AZIMUTH 

132 deg 

4.4 mdeg 

1 . 5 mdeg 

66 

ELEVATION 

15 deg 

6.7 mdeg 

2 . 7 mdeg 

60 

RANGE 

27 kft 

5.5 ft 

1.7 ft 

69 

LOKI SPHERE 
(RADAR #5) 
(ORIGINAL DATA) 
AZIMUTH 

137 deg 

2.9 mdeg 

2.5 mdeg 

14 

ELEVATION 

42 deg 

4.7 mdeg 

3.8 mdeg 

19 

RANGE 

27 kft 

6.0 ft 

3.5 ft 

42 

LOKI OPTICAL 
(RADAR # 3) 
(ORIGINAL DATA) 
AZIMUTH 

139 deg 

11.4 mdeg 

9 . 5 mdeg 

17 

ELEVATION 

78 deg 

11.4 mdeg 

8.2 mdeg 

28 

RANGE 

220 kft 

6.8 ft 

2.0 ft 

71 


in elevation. This difference may in part be caused by the 
much larger nominal slant range for the Optical payload as 
compared to the Zuni. The rms accuracy for Radar No. 5, 
tracking a Loki Sphere, is similar to the results for Radar 
No. 3 tracking the Zuni trajectories. 

3.5.2 Rms Errors in Positional Data Products 

The rms error estimates for positional data products 
are presented in Table 3.5-2. The rms latitude errors are in 
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TABLE 3.5-2 


ESTIMATED RMS OF NOISE-LIKE ERRORS IN 
POSITIONAL DATA 


DATA SET 

RMS 

ERROR 

MEAN 

SLANT 

RANGE 

ZUNI 

(RADAR It 3) 
LATITUDE 
LONGITUDE 
HEIGHT 

3 pdeg 
5 pdeg 
0.4m 

27 kft 

LOKI SPHERE 
(RADAR #5) 
(ORIGINAL DATA) 
LATITUDE 
LONGITUDE 
HEIGHT 

6 pdeg 
8 pdeg 
0.8 m 

27 kft 

LOKI OPTICAL 
(RADAR It 3) 
(ORIGINAL DATA) 
LATITUDE 
LONGITUDE 
HEIGHT 

65 pdeg 
73 pdeg 
2.1 m 

220 kft 


the range from 3 pdeg to 65 pdeg, with the largest error oc- 
curring at the largest slant range. The rms longitude errors 
are similar and range from 5 pdeg to 73 pdeg. Estimated rms 
height errors range from 0.4 m (1.3 ft) to 2.1 m (6.9 ft). 
Again, the largest error occurs at the largest slant range. 


3.6 SMOOTHING-FILTER TRAJECTORY BIAS 

The smoothing filters are intended to reduce high- 
frequency noise levels in the radar tracking data without 
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significantly changing the low-frequency signal component 
representing the payload motion along the trajectory. In prac- 
tice, the smoothing filter systematically distorts the nominal 
trajectory slightly. This section describes the results of an 
analysis of this systematic error component in smoothed tracking 
data . 

The technical approach of this analysis is explained 
with the aid of Fig. 3.6-1. The raw tracking data (azimuth, 
elevation, or range measurements) are decomposed into two signal 
components: a polynomial that represents the nominal trajectory 

and the residual data that represent the noise-like errors. 
This decomposition was previously introduced in Chapter 2. 

The effect of the smoothing filter is to smooth the residual 
data and to distort slightly the nominal trajectory polynomial. 
The distortion of the polynomial is termed smoothing- filter 
trajectory bias because it can produce a systematic bias-like 
error in the smoothed data. 


R 98660 


RAW TRACKING DATA 



SMOOTHED TRACKING DATA 

NOMINAL 

TRAJECTORY 

POLYNOMIAL^"" 



DISTORTED 

NOMINAL 

TRAJECTORY 

POLYNOMIAL 





WFF 

SMOOTHING 

FILTER 



+ 



+ 

"ESIDUAL 



SMOOTHED 

... 



RESIDUAL "*'** 
DATA 



Figure 3.6-1 Smoothing-Filter Trajectory Bias 
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An example of smoothing- filter trajectory bias is 
shown in Fig. 3.6-2. The plot on the left depicts the raw 
range data from Radar No. 3 for a Super Loki Optical payload. 
To model the trajectory signal, an orthogonal polynomial was 
fitted to the range data using the technique discussed in 
Section 2.2. The plot in Fig. 3.6-2 labeled "Range-Bias" shows 
the distortion of this polynomial that was produced by passing 
it through the WFF F00.40.10 smoothing filter. More than 
3 ft of bias was produced over most of the data segment. This 
bias is larger than the estimated rms noise level (2.0 ft) in 
the smoothed data. 

Trajectory bias from the smoothing filter does not 
always produce errors larger than the noise. Figure 3.6-3 
shows the azimuth and elevation biases for the Loki Optical 
data from Radar No. 3. These biases are much less than the 
rms estimated noise levels (9 mdeg). 


SLANT RANGE 


RANGE BIAS 



Filter-Induced Range Bias for Loki Optical 
Tracking Data, Radar No. 3 
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Figure 3.6-2 





AZIMUTH (mdeg) 


AZIMUTH 8IAS 



ELEVATION BIAS 

ft-08658 



TIME (sec) 


Figure 3.6-3 Filter-Induced Azimuth and Elevation Bias 

for Loki Optical Tracking Data, Radar No. 3 

The analysis of trajectory bias leads to the conclu- 
sion that the smoothing filters can produce systematic errors 
that are larger than the rms noise levels in the smoothed data . 
A recommended way of avoiding this error is to smooth only 
residual tracking data and then add the trajectory polynomial 
to the smoothed residuals at the output of the filter. 
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3.7 


SUMMARY AND CONCLUSIONS 


3.7.1 Summary of Task 2 

Under Task 2 the following objectives were met: 


• The propagation of noise- like errors in 
radar tracking data into positional data 
products was analyzed. A state-space 
covariance analysis was performed based 
on the stochastic error models developed 
under Task 1 . 

• The systematic smoothing error, termed 
smoothing- filter trajectory bias, was 
identified and analyzed. 


3.7.2 Conclusions 


The analysis results of Task 2 lead to the following 
main conclusions: 




Rms noise-like errors in smoothed track- 
ing data vary with payload for radars 
No. 3 and 5. The estimated rms noise 
levels of the smoothed data and positional 
data products are in the following ranges 
for the data analyzed in this study: 


AZ and EL : 
RANGE : 

LAT and LONG: 


1.5 mdeg to 9.5 mdeg 
1.7 ft to 3.5 ft 
3 Mdeg to 73 Mdeg 


HEIGHT: 1.3 ft to 6.9 ft 


• The bias-like errors caused by WFF smooth- 
ing filters can exceed the rms noise 
levels in smoothed tracking data when 
raw tracking data are passed through the 
filters . 
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TASK 3 - ALTERNATIVE NOISE REDUCTION TECHNIQUES 


4 . 1 INTRODUCTION 

The purpose of this chapter is to present recommenda - 
tions for improving the noise reduction techniques current! 
used at WFF for processing radar tracking data . These recom- 
mendations are based on the results of Tasks 1 and 2 of this 
study and an assessment of practical alternative algorithms. 

The current WFF data processing system is depicted in 
Fig. 4.1-1. At the top of this figure, the radar data tape 
provides inputs to the PASS 1 program, which produces as its 
output a working data tape. Data calibration and editing (e.g., 


RADAR 
DATA TAPE 



POSITIONAL 
DATA TAPE 


PASS 1 PROGRAM 


DATA PROC. REAP, ZONBtT PROGRAMS 


SMAD PROGRAM 


MESUP AND POSDAT PROGRAMS 


Figure 4.1-1 Review of Current WFF Smoothing Technique 
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to correct extreme data outliers) are then performed using 
programs such as DATA PROC and ZONBIT. The calibrated and 
edited data are next processed using the SMAD program to reduce 
the noise level in the data. It is in the SMAD program that 
the smoothing filters (analyzed under Task 2 of this study) 
are used to process the azimuth, elevation, and range data. 
The smoothed data are finally processed using the MESUP or 
POSDAT programs to produce data products, including data on 
payload position expressed in latitude, longitude, and height. 
The rms noise levels in these positional data products were 
estimated in Task 2 of this study. 


A. 2 DISCUSSION OF CURRENT FILTERING TECHNIQUES 

The current WFF smoothing filters are low-pass zero- 
phase finite-impulse-response (FIR) filters. The parameters 
of these filters are specified each time the SMAD program is 
run. Appropriate parameter values depend on the radar, the 
type of data being processed (azimuth, elevation, or range), 
the payload, and the portion of the trajectory that is being 
estimated. For example, data for ECC BALLOON, SUPER LOKI 
OPTICAL, and ZUNI ROCKETS are usually smoothed with filter 
parameters designated by the WFF code "F00 . 040 . 10 . " In contrast, 
THRUSH data may be analyzed using an F00.090.04 filter. For 
SCOUT and TAURUS ORION data, a variety of different filter 
parameters may be used, corresponding to different segments of 
the trajectory. 

According to the error analysis presented in Chapter 3, 
the rms noise levels in ZUNI, SUPER LOKI SPHERE, and SUPER LOKI 
OPTICAL data are significantly reduced by filter F00.040.10. 
However, the percent reduction of the noise level varied from 
14% to 71% for the data sets analyzed in this study. A main 
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conclusion to be reached from these results is that the current 
VFF filters can be effective in reducing rms noise levels, but 
the final noise levels in the smoothed data vary from data set 
to data set, and these noise levels are not estimated by the 
current VFF smoothing techniques . 

Another finding from Chapter 3, is that the current 
VFF smoothing technique can result in a systematic error, 
called filter-induced trajectory bias . This bias-like error 
is a slight distortion of the nominal trajectory signal as it 
is processed by the smoothing filter . For some of the data 
analyzed, the bias exceeded the rms noise level of the smoothed 
radar data. This finding can be interpreted positively as a 
verification that the noise levels are currently low and do 
not have to be reduced. Or the bias can be viewed as a known 
error source that should be eliminated. Fortunately, this 
error can be avoided by a simple modification of the current 
WFF smoothing procedure, as explained in Section 4.3. 

4.3 ALTERNATIVE NOISE REDUCTION TECHNIQUES 

4.3.1 Avoiding Filter-Induced Trajectory Bias 

The smoothing filter causes trajectory bias because 
in the SMAD program the nominal trajectory signal is smoothed 
along. with the noise in the data. The nominal trajectory signal 
is much larger than the noise, and so very small relative dis- 
tortions of the trajectory signal cause bias-like errors in the 
filter output. To avoid this bias, a polynomial estimate of the 
nominal trajectory signal should be subtracted from the radar 
data to yield residual radar data. The residual data, which 
consist mostly of measurement noise and only small-scale payload 
motions, would be smoothed using the current WFF filtering 
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technique. Finally, the polynomial estimate of the nominal tra- 
jectory signal would be added back into the smoothed residual 
data. By using this technique, the current WFF filters can 
reduce the noise level without producing significant bias error. 

4.3.2 Estimating and Minimizing Noise Levels in 
Smoothed Data 

As discussed in Section 4.2, the current WFF smoothing 
technique yields no estimate of the noise levels in the smoothed 
radar data. As a consequence, the noise levels in positional 
data products are also left unestimated. To correct these 
deficiencies, a more complicated data processing algorithm is 
required. In this section, two alternative approaches to esti- 
mating noise levels are discussed. The first is the simpler, 

.and requires an additional stage of signal processing to estimate 
the high-frequency noise levels in the data. The second approach 
is capable of higher accuracy, but is much more complicated. 

It is based on Kalman smoothing techniques, which are optimal 
with respect to prior information about the geometry of the 
tracking system, the physics of the descending payload, and 
the statistics of the measurement noise. 

Estimating High-Frequency Noise Levels - The high- 
frequency noise levels in smoothed radar tracking data (out- 
puts from the current SMAD program) can be estimated using the 
following procedure: 

• Step 1 - Select segments of the tracking 
data for analysis . The lengths of these 
segments may be as short as a few hundred 
data samples (e.g., ZUNI data processed 
under Task 1 of this study typically 
contained about 300 to 400 measurements 
per trajectory). For long trajectories 
spanning tens of thousands of measurements, 
segments may be selected from the beginning, 
middle, and end of the trajectory. 
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• Step 2 - For each segment compute residual 
azimuth, elevation, and range data by sub- 
tracting orthogonal polynomials from each 
channel of data . The polynomials are 
fitted to each data segment using least- 
squares as discussed in Appendix A. The 
degree of each polynomial is not critical 
for estimating high-frequency noise levels. 
Appropriate polynomial degrees for the 
tracking data analyzed under Tasks 1 

and 2 of this study are in the range 
6 to 10 for segments containing about 
500 data samples. 

• Step 3 - Process each channel of residual 
data with the autoregressive modeling 
algorithm discussed in Appendix B . The 
output of this procedure is a stochastic 
model for the residual data. 

• Step 4 - Use the autoregressive models 
from Step 3 to estimate the power spectra 
of the azimuth, elevation, and range data 
for each segment" ! The rms high-frequency 
noise level in each data set is inferred 
from the level of its power spectrum at 
high-frequencies as indicated in Fig. 4.3-1. 


This method for estimating noise levels can be implemented in 
a new computer program. By running this program on the outputs 
of SMAD, POSDAT, or MESUP, noise levels can be estimated for 
smoothed tracking data or positional data products. 

An Alternative Smoothing Technique - Optimal estimation 
techniques can be used instead of the current WFF smoothing 
procedure. The motive for using optimal smoothing is to obtain 
the most accurate data products together with reliable error 
estimates. This requires that the smoothing algorithm be opti- 
mized with respect to several kinds of prior information about 
the tracking system geometry, the physics of the descending 
payload, and the statistics of the measurement noise and other 
uncertainties. Incorporating this information optimally requires 



A 7556 


ESTIMATED 



Figure 4.3-1 Estimation of High-Frequency Noise Variance 

from Power Spectral Density (PSD), 
f = Data Sampling Frequency (Hz) and 

S Q = Nominal High-Frequency Noise PSD 

a much more complicated processing algorithm than the one cur- 
rently used at WFF. In the majority of cases, where well cal- 
ibrated low-noise tracking data are being processed, the added 
complexity of an optimal estimator may not be justified. How- 
ever, for special cases in which the measurement noise is 
large, the radar calibration is not precise, or the tracking 
data contain long gaps due to missing measurements, optimal 
smoothing may be essential for meeting the objectives of the 
data analysis. 

The recommended way of implementing an optimal esti- 
mator for processing radar tracking data is to use Kalman 
filtering and smoothing algorithms. There is a well developed 
software technology and mathematical theory to support the de- 
sign and implementation of Kalman-type processors (e.g. , Refs. 5 
and 19). Moreover, the mathematical formalism is very . flexible , 
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which makes it possible to improve the accuracy of the processor 
by adding new or improved error models or additional data inputs 
to the smoother, without redesigning the software. 

Specific examples of tracking data that are candidates 
for the use of Kalman estimation techniques are NIKE-ORION and 
TERRI ER-MALEMUTE data from Radars Nos. 8 and 41 in Peru. For 
these examples, the payloads were tracked above the atmosphere 
(above 50 km) so that the physics of the payload descent can 
be represented with a relatively simple mathematical model. 

The primary sources of uncertainty in the estimates of payload 
position are uncertainties in the radar calibration, random 
noise in the radar data, and possible data gaps caused by missing 
measurements or high noise levels. A Kalman smoothing algorithm 
for these data sets is developed and verified under Task 4 of 
this study. The detailed specification of the algorithm and 
examples of its performance with data provided by WFF is pre- 
sented in Chapter 5. 

4.4 SUMMARY OF RECOMMENDATIONS FOR IMPROVING THE CURRENT 
WFF NOISE REDUCTION TECHNIQUES 

Based on the results of (1) the error analyses of the 
current WFF radar data smoothing fiters (conducted under Tasks 1 
and 2 of this study) and (2) an assessment, of alternative tech- 
niques, the following recommendations are made for improving 
the noise reduction techniques used at WFF: 

• To avoid filter-induced trajectory bias , 
residual radar data should be smoothed 
instead of raw radar data. The residual 
data may be computed by subtracting from 
the raw data low-degree orthogonal poly- 
nomials, which are least- squares estimates 
of the nominal trajectory signal in the 
raw data. 
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• To estimate the high-frequency rms noise 
levels in smoothed tracking data produced 
using the existing smoothing filters , 
residual smoothed data ( i . e . , either 
smoothed residual data or smoothed raw 
data minus an estimated nominal trajec- 
tory signal) may be processed with the 
autoregressive (AR) modeling algorithm 
discussed in Appendix B. The noise level 
is estimated from the high-frequency 
part of the power spectrum of the AR 
model . 

• To process the radar tracking data opti- 
mally, an alternative Kalman filter/ - 
smoother algorithm is recommended . Optimal 
estimation is much more complicated than 
the current smoothing procedure because it 
uses prior information about tracking 
system geometry, the physics of the de- 
scending payload, and the statistics of 
the measurement noise, radar calibration 
errors, and other uncertainties. It is 
expected that optimal • smoothing tech- 
niques would be most appropriate for 
processing tracking data having high 
noise levels, significant data gaps, or 
large uncertainties on radar calibration 
errors. 
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5 . . TASK 4 - DEVELOPMENT AND VERIFICATION OF AN 

ALTERNATIVE SMOOTHING ALGORITHM 


This chapter describes the development and verifica- 
tion of a new smoothing algorithm for processing radar tracking 
data . The algorithm takes as inputs the tracking data (azimuth, 
elevation, and range) from one radar and provides as outputs 
estimates of payload position and velocity as functions of 
time. The algorithm is a Kalman filter/ smoother that is optimal 
(i.e., unbiased and minimum-variance) with respect to prior 
information about the tracking system geometry, the physics of 
the payload descent, and the statistics of measurement noise 
and radar calibration errors. The algorithm, in its present 
form, is intended for estimating trajectories above the atmos- 
phere (height > 50 km). The algorithm detects and appropriately 
processes isolated data outliers. Moreover, extended measure- 
ment gaps (caused by missing data or high noise levels) are 
processed optimally when their locations in the data set are 
specified as input parameters. The performance of the algorithm 
is verified using tracking data provided by WFF for this 
investigation. 

This chapter is organized- as follows. Section 5.1 
describes the Task 4 objectives and the technical approach for 
meeting these objectives. Section 5.2 provides an overview of 
the algorithm, while the mathematical details are discussed in 
Sections 5.3 through 5.7. The verification of algorithm per- 
formance with tracking data from Radars Nos. 8 and 41 in Peru 
is presented in Section 5.8. The accomplishments of this task 
are summarized in Section 5.9. 
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5.1 


INTRODUCTION 


5.1.1 Obj ective 


Task 4 has two main objectives: 


• Develop an algorithm for processing noisy 
radar tracking data to estimate payload 
trajectories above 50 km in altitude. 
The algorithm should handle data outliers 
and gaps caused by missing data or high 
noise levels, and should provide a real- 
istic estimate of rms error for the esti- 
mated trajectory. 

• Verify the performance of the algorithm 
using radar data provided by WFF. 


5.1.2 Technical Approach 


To meet the two objectives, the algorithm was devel- 
oped using the established theory of Kalman optimal filtering 
and smoothing (Ref. 5). The technical approach consists of 
four steps: 


• Compute a Nominal Trajectory - The nom- 
inal trajectory is computed using (1) a 
best-guess initial position and velocity 
for the payload, and (2) deterministic 
models for normal gravitation and nominal 
atmospheric drag accelerations'. 

• Compute Nominal Radar Measurements - The 
nominal radar measurements are time series 
for azimuth, elevation, and range corres- 
ponding to an ideal radar tracking the 
nominal trajectory. 

• Compute Residual Radar Measurements - The 
residual radar measurements are defined 
as the actual measurements minus the 
nominal measurements computed in Step 2. 
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• Estimate Corrections to the Nominal 

Trajectory - The corrections to the nom- 
inal trajectory are computed by processing 
the residual radar data from Step 3 using 
a Kalman f il ter/ smoother . The final opti- 
mal trajectory estimate is then computed 
by adding these corrections to the nominal 
trajectory from Step 1. The error covar- 
iance matrices of the trajectory estimates 
are computed by the Kalman smoother. 

There are two main advantages to processing residual 
radar data. The first is that the relation between residual 
tracking data and residual payload motions about a nominal 
trajectory can be accurately modeled by linear time-varying 
state-space difference equations . (The reason for this is 
that, with high probability, the actual payload motion is a 
small percent perturbation about the nominal trajectory. This 
expectation is justified for trajectories above the atmosphere, 
i.e., height > 50 km, because at high altitudes the accelera- 
tions caused by atmospheric drag and gravitation, can be ade- 
quately represented using simple models.) The state-space 
difference equations are precisely the type of mathematical 
model that is consistent with the recursive Kalman filter/ 
smoother algorithms. 

The second advantage to processing residual radar 
data is that they can be accurately modeled as realizations 
of zero-mean non-stationary random processes . Moreover, these 
processes can be represented as outputs from linear time-varying 
state-space difference equations driven by white noise. This 
is the type of data for which Kalman filter/smoother algorithms 
are statistically optimal (i.e., the trajectory estimates are 
unbiased and have the smallest possible error variances). 
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5.2 


OVERVIEW OF ALGORITHM 


The front end of the trajectory estimation algorithm 
is a pre-processor that (1) computes a nominal trajectory for 
the payload, and (2) processes the raw radar data (e.g., output 
data from the WFF PASS 1 program) to produce residual radar 
data as outputs. Figure 5.2-1 depicts the pre-processor, its 
inputs, and its outputs. The inputs are: 


• Best-Guess Initial State of Payload - The 
estimated position and velocity of the 
payload at the initial time expressed in 
earth-centered Cartesian inertial 
coordinates 

• Nominal Gravitation Model - The normal 
corrections to a point-mass gravitation 
model, which account for the oblateness 
of the earth's gravitational field 

• Nominal Atmospheric Drag Acceleration 
Model - The expected atmospheric drag as 
a function of altitude and payload veloc- 
ity, based on the U.S. Standard Atmosphere, 
1976 (Ref. 6) 

• Geodetic Coordinates of the Tracking Radar - 
Expressed in terms of latitude, longitude, 
and height with respect to the reference 
ellipsoid currently used in the WFF data 
processing programs 

• Actual Radar Measurements - Time series 
of azimuth, elevation, and range measure- 
ments taken at uniformly spaced time 
intervals . 


The outputs of the pre-processor are used by the 
filtering and smoothing stages of the algorithm. As indicated 
in Fig. 5.2-1, the outputs are: 



Figure 5.2-1 Block Diagram of Pre-Processing for 

Trajectory Estimation 


• Keplerian Payload Trajectory - The ideal- 
ized trajectory which the payload would 
follow if the earth were a point mass, 
the atmospheric drag accelerations were 
zero, and the initial position and veloc- 
ity of the payload were known exactly 

• Nominal Payload Trajectory - The expected 
trajectory for the payload , given the 
initial estimate of the payload position 
and velocity at the initial time, the 
normal gravitation of the earth, and a 
model for the nominal drag accelerations 
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• Residual Radar Measurements - The trans- 
formed radar data, which contain all 
available information about the actual 
departure of the payload trajectory away 
from the nominal trajectory. The residual 
data are inputs to the Kalman filter in 
the second part of the trajectory esti- 
mation algorithm. 


The second half of the trajectory estimation algorithm 
is a Kalman filter/smoother. A block diagram of the filter 
and smoother is presented in Fig. 5.2-2. As shown in this 
diagram, the Kalman filter algorithm has six inputs: 


• The Processing Mode - a parameter that 
determines whether the algorithm is to 
be optimal for base-line measurement 
noise (mode 1), for data gaps caused by 
missing data or very noisy data (mode 2), 
or for automatic detection and optimal 
processing of isolated data outliers 
(mode 3) 

• An Estimate of the Payload’s Initial State 
Vector and the Error Covariance Matrix of 
this Estimate - the best-guess estimate 
(••before any radar data are processed) of 
the payload's position and velocity with 
respect to the nominal trajectory at the 
initial time 

• The Residual Radar Measurements and the 
Time Between Successive Measurements - the 
residual radar data computed by the pre- 
processor 

• The Keplerian Trajectory - the trajectory 
based on a point-mass earth, which is 
computed by the pre-processor 

• The Nominal Radar Measurements - nominal 
tracking data corresponding to an ideal 
radar tracking the nominal payload tra- 
jectory, which are computed by the pre- 
processor 
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• The Error Model - a state-space stochastic 
model for radar measurement noise and 
radar calibration errors. (The model 
could also represent acceleration noise 
caused by small unpredictable perturba- 
tions in the gravitational field away 
from the normal field for an ellipsoidal 
earth. Based on the results of processing 
WFF radar data from Peru for Radars Nos. 8 
and 41, it was concluded that gravitational 
errors are much smaller than radar measure- 
ment errors. Therefore, an error model 
for gravitational noise is not included 
in the present version of the trajectory 
estimation algorithm.) 


The residual radar measurements are processed causally 
by the Kalman filter, starting with the initial time and pro- 
ceeding to the end of the data set. The outputs of the Kalman 
filter are three time series: 
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• Filtered State Estimates - The estimated 
position and velocity of the payload 
with respect to the nominal trajectory 
as a function of time, together with 
additional state variables that model 
radar measurement errors. These esti- 
mates are optimal with respect to a causal 
processing of the data. 

• Error Covariances - The error covariance 
matrices of the filtered state estimates. 

• Innovations and Their Variances - The 
innovations are the prediction errors 
made by the filter as it predicts what 
the next residual measurements of azimuth, 
elevation, and range will be one time-step 
ahead. In the data-adaptive mode (mode 3) 
of the trajectory estimation algorithm, 
the innovations and their variances are 
used to detect outliers. 


As indicated in Fig. 5.2-2, the smoothing algorithm 
processes the outputs of the Kalman filter. The smoother is 
also recursive, but it works backward in time, starting at the 
end of the input time series, and running back to the initial 
time. Because the smoother runs in reverse, the outputs of 
the filter are stored in random access files so that they can 
be accessed in the reverse order to which they were stored. 


The outputs of the smoother are optimal estimates of 
the payload position and velocity (with respect to the nominal 
trajectory) as a function of time and the error covariances 
for these estimates. Also included in the outputs of the smoother 
are estimates of the radar calibration errors that are modeled 
by the algorithm and their covariances. 


In Sections 5.3 through 5.7, the mathematical details 
of the trajectory estimation algorithm are discussed. The 
organization of these sections follows the order in which cal- 
culations are performed when radar data are being processed. 
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5.3 KEPLERI AN -TRAJECTORY MODULE 

The first step in computing the nominal payload tra- 
jectory is to compute a good approximation to it, called the 
Keplerian trajectory . The Keplerian trajectory is the position 
and velocity of the payload as a function of time for the ideal- 
ized model of a point-mass earth and a point-mass payload. 
Using earth-centered Cartesian inertial coordinates (defined 
in Appendix D), the payload position r(t) and velocity v(t) at 
time t can be expressed as linear combinations of the position 
and velocity at t=0. 


v(t) = r(t) = velocity vector 

(5.3-1) 

r(t) = f ( t ) r ( 0 ) 

+ g(t)v(0) 

(5.3-2) 

v(t) = fe( t )r(0 ) 

+ g(t)v(0) 

(5.3-3) 


In Eqs. 5.3-2 and 5.3-3, the scalar functions f(t) and g(t) are 

known as the " f and g functions" (Refs. 7 and 20). They are 

computed using the canonical units of length and time, UL and UT, 

which are defined in terms of the semi-major axis, a (m] , of 

the reference ellipsoid for the earth and the gravitational 
3-2 

constant GM (m *s ] of the earth as follows: 


UL = a 



The initial position r(0) and 
the canonical units: 


[m] (5.3-4) 

l s ) (5.3-5) 

velocity v(0) are scaled using 
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r Q = r(0)/UL 

(3.5-6) 

v = v(0) *UT/UL 
— o — 

(3.5-7) 


From these scaled initial conditions, the following parameters 
are computed (superscript T denotes a matrix transpose): 


l£ol = 



(5.3-8) 

IX 0 ! 2 = 

T 

V V 

— o — o 


(5.3-9) 

RV = 

T 

r v - 
— o — o 

2 

(5.3-10) 

a = 

l< 

o 

ro 

i 


(5.3-11) 

c o = 

1 + a • | 

So' 

(5.3-12) 

A = 

<-a>^ 


(5.3-13) 


The parameters defined by Eqs. 5.3-8 to 5.3-13 are used to 
compute the four additional parameters E q , e, M q , and N as 
follows: 

If c / 0 then E = tan" 1 l-^L] + n-(c < 0) (5.3-14) 

o o c*A o 

o 

where the quantity n*( c Q < 0) = n if c q < 0, and 0 otherwise. 

If c Q = 0 then E q = sign(RV)*n/2 (5.3-15) 

where sign(RV) = 1 if RV > 0 
sign(RV) = 0 if RV = 0 
sign(RV) = -1 if RV < 0 
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(5.3-16) 


If cos(E q ) / 0 then e = c q /cos(E o ) 


If cos(E q ) = 0 then e = 


A* sin(E Q ) 


M = E - RV/A 
o o 


N = A -3 


For each time t of interest, the f and g functions 
and their derivatives are computed as follows: 


f(t) = 1 - (1 - cos ( E - E o ])* t ^ t 

1 — o 1 


g(t) = [t/UT - IE - E q - sin(E - E Q ) J/n] 


t(t) = ['7TFT • sln( 


E - E. 


1/UT 


g(t) = 1 - ~ • [1 - cos ( E - E 0 )] (5.3-23) 

In Eqs. 5.3-20 to 5.3-23, the parameters r and E are computed 
using the following three-step iterative technique: 


Set M = M + N- t/UT 
o 


£ = 10 ' 
max 


E = 0 


(5.3-26) 


Step 1 set e = M + e*sin(E) - E 


(5.3-27) 
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(5.3-28) 


(5.3-29) 

The algorithm specified by expressions 5.3-24 to 5.3-29 is the 
recommended way of solving Kepler's equation for E: 

M - E - e-sin(E) (5.3-30) 

5.4 NOMINAL -TRAJECTORY MODULE 

This section discusses the mathematical details of 
computing a nominal trajectory for the payload at altitudes 
greater than 50 km. The nominal trajectory is computed by 
adding small corrections to the Keplerian trajectory defined 
in Section 5.3. These corrections are based on deterministic 
models for the normal gravitation of the earth and the expected 
atmospheric drag force acting on the payload. 

The approach is to compute the perturbations in gravi- 
tation and the small atmospheric drag forces which the payload 
would experience if it were to follow the Keplerian trajectory. 
These small quantities are then used to compute the corrected 
payload trajectory (nominal trajectory) by using the linearized 
equations of motion for small departures from the Keplerian 
trajectory. 

5.4.1 Normal Gravitation 

The normal gravitation of an ellipsoidal earth model 
is computed using the reference ellipsoid parameters listed in 


Step 2 if I e I > e 
c — 1 ' max 

then set E = E + e/1.4 
and go to Step 1 

else , go to Step 3 

Step 3 set r = A 2 *[l - e*cos(E)) 
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Table 5.4-1. These are the parameter values for the reference 
ellipsoid used in the WFF POSDAT program. The normal gravita- 
tional acceleration a (r) at position r with respect to the 
center of the reference ellipsoid is expressed as follows: 

a (r) = - G — o • r + 6a (r) (m*s~ 2 ] (5.4-1) 

g |r| J g ' 

The first term in Eq. 5.4-1 is the gravitation of a point-mass 
earth model. The second term is the correction for the normal 
gravitation of the reference ellipsoidal model. The radial 
and tangential components of the correction are first computed, 
then they are transformed to Cartesian inertial coordinates.' 
The equations for these calculations are given below: 

6a radial = <3/4)C 2 G 2 [l + 3-cos(2*A 0 )] [m-s' 2 ] 

(5.4-2) 

6a tangnt = ~ ^ 3/4 )C 2 G 2 [ c ° s < 3 • a q ) - cos(A q ) ]/ sin(A Q ) 

(5.4-3) 


TABLE 5.4-1 

GEODETIC PARAMETERS FOR REFERENCE ELLIPSOID 


a = 

6.378166* 10 6 [m] 

= semi-major axis 

b = 

6.356784-10 6 [m] 

= semi-minor axis 

f = 

1/298.3 

= flattening 

Q = 

7 . 292115147 • 10~ 5 [rad/s] 

= earth's angular 
rotation speed 

GM = 

3.986005-10 14 [m 3 *s" 2 J 

= earth's gravitational 
constant 
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In Eqs. 5.4-2 and 5.4-3, the parameters C 2 > G, and A q are com- 
puted as follows: 

C 2 = 1.082626* 10~ 3 /(GM/a 2 ) (5.4-4) 

G = GM/ I r | 2 (5.4-5) 

A q = (n/2) - A ^ (5.4-6) 

The parameter A^ (north latitude) in Eq. 5.4-6 is computed 
using the earth-centered Cartesian inertial coordinates r of 
the payload position: 


r = 


, i 

r 2 r 3 

(5.4-7) 

r l2 = 


+ 4 

(5.4-8) 

III = 

$ 

+ 4 + r 3 

(5.4-9) 

If r^ 2 '> 0 then A^ = 

tan 

1 ( r 3 / r 12 ) 

(5.4-10) 


If r ^2 = 0 then A^ = sign( r^ ) • n/2 (5. 4-11) 

In Eq. 5.4-11, the sign(x) function is defined as follows: 
sign(x) = 1 if x > 0; sigri(x) = -1 if x < 0; and sign(x) = 0 
if x = 0. 

The Cartesian inertial components of the gravitational 
correction are computed using the radial and tangential compo- 
nents (defined in Eqs. 5.4-2 and 5.4-3) as follows: 

6a = [6a , 6a 0 6a 0 ]^ (5.4-12) 

-g gl g2 g3 

6a g i = G 12 *CL (5.4-13) 
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G, ~ * SL 


(5.4-14) 



In Eqs. 5.4-13 to 5.4-15, the parameters G-^ » CL, and SL, are 
defined as follows: 


°12 = 6a radial' cos(A l > + 6a tangnt‘ sin<A l> (5.4-16) 

If r ^2 > 0 then CL = r^/r^ (5.4-17) 

and SL = * 2^12 (5.4-18) 

If r 12 = 0 then CL = 0 (5.4-19) 

.and SL = 0 (5.4-20) 

5.4.2 Atmospheric Drag 


The nominal atmospheric drag force f^ acting on the 
payload at altitudes greater than 50 km is computed as follows: 

il< h, -p/a^ = ' 0/2 ) * c d ’ (h) * 1 Vp/a 1 • y p /a (5.4-21) 

The symbols in Eq. 5.4-21 have the following meanings: 


id 

- 

drag force (vector) 

[N] 

h 


payload altitude 

lm] 

-p/a 

= 

velocity of payload with 
to the atmosphere 

respect , 

[m* s~ 

C d 

= 

payload drag coefficient 

[1 < C 

A P 

= 

2 

payload cross-sectional area [m ] 

P (h) 

= 

atmospheric mass density 

[kg*m" 
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The nominal mass density is computed using the follow- 
ing model, which is based on the U.S. Standard Atmosphere, 
1976 (Ref. 6): 

p (h ) = 2.2* 10’ 6 ‘ 5 * 10 5 * h (5.4-22) 


In Eq. 5.4-22, the payload height h above the reference 
ellipsoid is computed using the ellipsoid parameters (a = semi- 
major axis, f = flattening, b = (l-f)*a = semi-minor axis) and 
|r| , the distance of the payload from the center of the ellipsoid 


h = [r 12 ]/cos (LT) ] - N 

LT = . tan 1 (T/(l - f) 2 ] 

T = r 3 /|r| 

N = a 2 - (a 2 -cos 2 (LT) + b 2 - sin 2 ( LT) ] ~ 1/2 


(5.4-23) 

(5.4-24) 

(5.4-25) 

(5.4-26) 


In Eq. 5.4-21, the nominal velocity Yp/ a °f the payload 
with respect to the atmosphere is computed from the position 
of the payload r with respect to the center of the earth, the 
velocity v of the payload with respect to an inertial frame, 
and the angular velocity 0 of the earth with respect to an 
inertial frame: 

Vp/ a = v - Qxr (5.4-27) 


The three Cartesian inertial components of Yp/ a * an< * 9- are 
indicated in the following equations: 


V . = [v'V v (2 > v< 3 >) T 

-p/a p/a p/a p/a 


(5.4-28) 


v 



v. 



(5.4-29) 


Q = [0 0 «] T 


(5.4-30) 
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In Eq. 5.4-30, the first two components of ft are zero because 
the first two Cartesian inertial axes span the earth's equa- 
torial plane, while the third axis is the rotational axis of 
the earth model. (0 is the earth's rotational rate [rad*s ^].) 
Using Eqs. 5.4-28 to 5.4-30 in 5.4-27 leads to the following 
expressions for computing the relative velocity vector Y p / a : 


v 


( 1 ) _ 

p/a 


V 1 + 0 • r 2 


v 


( 2 ) _ 

p/a 


v 2 - n-r l 


V 


(3) _ 
p/a 


= V/ 


(5.4-31) 

(5.4-32) 

(5.4-33) 


By using Eqs. 5.4-22 to 5.4-33, with r and v evaluated 
on the Keplerian trajectory, the drag force vector f^ is computed 
from Eq. 5.4-21. The inertial acceleration a^ of the payload 
caused by the drag force f^ is then given as follows: 


^d 



(5.4-34) 


where m p is the payload mass (kg]. 

5.4.3 Trajectory Corrections 

In this section the equations are presented for com- 
puting a nominal payload trajectory , given the Keplerian tra- 
jectory defined in Section 5.3, the gravitational correction 
defined in Section 5.4.1, and the atmospheric drag acceleration 
defined in Section 5.4.2. The approach is to compute position 
and velocity corrections that are added to the Keplerian tra- 
jectory. These corrections are computed by solving the line- 
arized equations of motion for the payload, which govern small 
perturbations about the Keplerian trajectory. 
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The dynamical state of the payload on the nominal 
trajectory is represented by the state vector X, which is de- 
fined in terms of the nominal payload position R and velocity 
V. R and V are measured with respect to the center of the 
earth and are expressed in earth-centered Cartesian inertial 
coordinates : 

(5.4-35) 

R = IR 1 R 2 R 3 ) T (5.4-36) 

V = [V 1 V 2 V 3 ] T (5.4-37) 

The dynamical state of the payload along the Keplerian 
trajectory is defined in the same way: 




(5.4 -38) 


In Eq. 5.4-38, the position r and velocity v are expressed in 
earth-centered Cartesian inertial coordinates. The nominal 
and Keplerian state vectors are functions of time t and are 
related to each other as follows: 

X(t) = X (K) (t) + 6x(t) (5.4-39) 

In Eq. 5.4-39, 6x(t) is the correction that is added to the 
Keplerian state at time t to account for the influence of normal 
gravitation and atmospheric drag. In the remainder of this 
section, equations are presented for computing 6x(t) at uniformly 
spaced times t = k*6t, k = 0, 1, 2, ..., 6t = sampling inter- 
val [s]. The following notation is used for sampled values: 
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6x, - 6x(k*6t) for k = 0, 1, 2, ... 


(5.4-40) 


The linearized equations of motion for small perturba- 
tions about the Keplerian trajectory are as follows: 


6x(t) = F(t)6x(t) + 6u(t) 


(5.4-41) 


6u( t) 


6a(t) 


0 

= 0 
0 


(5.4-42) 


In Eq. 5.4-42, the disturbing acceleration 6a(t) is the sum of 
the gravitational and atmospheric drag accelerations defined 
in Sections 5.4.1 and 5.4.2: 


6a(t) = 6a g (t) +5a d (t) 

The matrix F(t) in Eq. 5.4-41 is defined as follows: 


O3 

F(t) = J 


(5.4-43) 


J(t) 0, 


(5.4-44) 


In Eq. 5.4-44, 0^ is the 3x3 zero matrix, I^ is the 3x3 identity 
matrix, and J(t) is the 3x3 matrix defined as follows in terms 
of the Keplerian position r(t): 


J(t) = 


|r(t)| 


r(t)lr(t)] 1 

3 , — - i. 

|r(t)r 


(5.4-45) 


In Eq. 5.4-45, GM is the gravitational constant multiplied by 
the mass of the earth, and r(t) is the position of the payload 
on the Keplerian trajectory. As shown in an unclassified sec- 
tion of Ref. 8, the sampled values of 6x(t), t = k*6t, k = 0, 
1, 2, ..., satisfy the following equation: 


6 -k+l " *k 6 -k + r k 6 -k 
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(5.4-46) 


In Eq. 5.4-46, the 6x6 transition matrix is defined as 
follows : 


= e 


. JF*6t _ *11 *12 


*21 *22 


(5.4-47) 


The 3x3 submatrices in Eq. 5.4-47 are computed using the follow- 
ing equations: 


|r(t)| 


(5.4-48) 


r( t ) • r ( t ) 
|r(t) | 2 


(5.4-49) 


*22 = l cosh ( a ‘^ t ’V2) - cos(a*6t)]-A + cosfa-fit)-!^ 


= 1/a * [ — ♦sinh(oc6t*>/2) 

. V2 


(5.4-50) 


sin(oc6t) ] -A + sin(a 


•6t)-l 3 


a • [ > /2*sinh(a*6t*V2) + sin(a*6t)] 


A - sin(a*6t) 


(5.4-51) 


(5.4-52) 


In Eq. 5.4-46, the 6x3 input weighting matrix is defined as 
follows : 


(5.4-53) 


^ = l/Q' 


cosh(a-6t* > /2)-l cos(a*6t)-l~| cos(cf*6t)-l 
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(5.4-54) 


y 2 =1 


sinh(ci’6t'-j2) _ sin(or6t) 
V2-a " 


sin(a’6t) 

A a x 3 


(5.4-55) 


The above equations are used to compute the state 
corrections 6x^ for k=0, 1, 2, given the disturbing 
acceleration 6a^. The sample values of the nominal tra- 
jectory are then computed by adding the corrections to the 
Keplerian state: 


X k = X< K) + «x k (5.4-56) 

'5.5 RADAR-MEASUREMENTS MODULE 
5.5.1 Introduction 


The purpose of this section is to present the mathemat- 
ical details of the radar-measurements software module. The 
module is used to compute both nominal radar measurements and 
residual radar measurements. As depicted in Fig. 5.5-1, this 
module has the following inputs: 

• Nominal Payload Trajectory - expected 
payload position and velocity expressed 
in earth-centered Cartesian inertial 
coordinates 

• Radar Position - geodetic coordinates 
(longitude , latitude, and altitude) of 
the tracking radar, defined with respect 
to the reference ellipsoid 

• Radar Tracking Data - real tracking data 
expressed in radar coordinates. 


The outputs of this module are (1) a time series of 
nominal radar measurements and (2) a time series of residual 
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Figure 5.5-1 Radar-Measurements Module 

radar measurements, both time series expressed in radar coor- 
dinates (azimuth, elevation, and range). The nominal measure- 
ments are the data that would be acquired if an ideal radar 
were to track the payload along the nominal trajectory. The 
residual measurements are the corrections that are added to the 
nominal measurements to produce the actual radar measurements. 


5.5.2 Transform to Geocentric Coordinates 


As indicated in Fig. 5.5-1, the inertial coordinates 
of the nominal trajectory and the geodetic coordinates of the 
radar are transformed to geocentric (earth-centered and earth- 
fixed) Cartesian coordinates. The equations governing these 
transformations are presented in the following: 
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The geodetic position coordinates of the tracking 
radar (LO = east longitude, LA = north latitude, H = height 
above reference ellipsoid) are transformed to geocentric 
Cartesian coordinates R^ ra< ^ ar ) using the following equations 
(a = semi-major axis and b = semi-minor axis of the reference 
ellipsoid) (Ref. 9): 
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(5.5-5) 


R (radar) _ 


R x = (N + H) •cos(LA)*cos(LO) (5.5-6) 

R 2 = (N + H) • cos(LA) • sin(LO) (5.5-7) 

R 3 = [(b/a) 2 *N + H]* sin (LA) (5.5-8) 

N = a 2 * [a 2 -cos 2 (LA) + b 2 -sin 2 (LA ) ]‘ 1/2 (5.5-9) 



The nominal payload position relative to the radar, 
j.(p/r) } i s expressed in geocentric Cartesian coordinates as 
follows: 


r (p/r) _ r (geo) _ R (radar ) 


(5.5-10) 


5.5.3 Transform to Topocentric Coordinates 

As depicted in Fig. 5.5-1, the next step in computing 
nominal radar measurements is to transform the payload position 
relative to the radar, to topocentric Cartesian coordinates 

R (ned) _ north, e = east, d = down) with the origin located 
at the radar position. The equations governing this transforma- 
tion are given in the following: 

p(ned) _ T (ned) (p/r) 

- ' T (geo) £ 

-sin(LA) •cos(LO) -sin(LA) • sin(LO) 

-sin(LO) cos(LO) 

-cos(LA) • cos(LO) -cos(LA) • sin(LO) 


T (ned) _ 
1 (geo) 


(5.5-11) 


cos(LA) 

0 

-sin (LA) 


(5.5-12) 
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R 


north 

(ned) _ 

n 

R e 

= 

east 


*d_ 


down 







•(5.5-13) 


5.5.4 Transformation to Radar Coordinates 


The final transformation indicated in Fig. 5.5-1 is 
the conversion of Cartesian north-east-down coordinates to 
spherical radar coordinates (AZ = azimuth, EL = elevation, 

RA = range). Azimuth is measured positive eastward, with 
AZ = 0 for due north. Elevation is measured positive toward 
the zenith, with EL = 0 for the horizontal. Range is measured 
positive away from the radar, with RA = 0 at the radar. The 
radar measurements are computed from R^ nec ^ using the following 
equations: 


R ne ‘ l- 

K + K 

n e 

(5.5-14) 

cos(AZ) 

= R /R 
n' ne 

(5.5-15) 

sin ( AZ) 

= V R ne 

(5.5-16) 

EL = tan 

' 1( -V R ne> 

(5.5-17) 

RA = Jr* 

2 2 ' 

+ R“ + Rj 
e d 

(5.5-18) 


5.5.5 Outputs 

The outputs of the radar-measurements module are the 
nominal measurement vectors (Z(t)) and the residual measurement 
vectors (z(t)) for each time t = k*6t, k = 0, 1, 2, ..., along 
the nominal trajectory: 
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Z k = Z(k* 6t ) 
z k = z ( k* 6 1 ) 



7 (actual) 7 
-k ' -k 


where 

^(actual) _ _ c . , . , . 

Z k = vector of actual tracking 

data at time t = k*6t 


(5.5-19) 

(5.5-20) 

(5.5-21) 


5.6 KALMAN -FILTER MODULE 
5.6.1 Introduction 

This section describes the mathematical details of 
the Kalman-filter module. The purpose of the filter is to 
process residual radar tracking data as inputs. The outputs 
are a time series of filtered estimates of the payload state 
(position and velocity) relative to the nominal trajectory 
defined in Section 5. A. Additional outputs are estimates of 
other state variables representing radar calibration errors, 
the error covariance matrices for the state estimates, and a 
time series of innovations data (the one-step-ahead prediction 
errors of the Kalman filter and their variances). 

The Kalman-filter algorithm presented in this section 
has several important properties (Ref. 5): 

• 

estimates of the payload state at time t. 

These estimates are based on (1) stochastic 
error models for the radar data, (2) the 
error covariance of the estimated initial 
state, and (3) the residual radar measure- 
ments up to, but not beyond, time t. 


The Kalman filter is a recursive algorithm 
for computing unbiased minimum-variance 


5-26 


• The filter is linear and time-varying 
because it estimates small corrections 
to the nominal dynamical state of the 
payload as a function of time along the 
nominal trajectory. 

• Missing radar data and isolated outliers 
are handled optimally . The filter proc- 
esses missing data by optimally extrapo- 
lating over the gaps (this is called 
mode 2 processing). In mode 3 processing, 
data outliers are detected automatically 
by comparing the innovations with their 
theoretical rms values. When the innova- 
tions exceed approximately three standard 
deviations, the noise variance of the 
measurement-noise model is continuously 
and automatically increased, and the new 
measurement is processed optimally with 
respect to this higher noise level. 

• The filter computes the theoretical rms 
accuracy of the estimated trajectory 
based on the radar error model, the uncer- 
tainty of the initial position and veloc- 
ity, and all radar measurements up to 
the current instant. 

• The outputs of the filter are sufficient 
statistics for computing the final 
smoothed estimate of the trajectory based 
on all available radar measurements. 


A block diagram of the data processing performed by 
the Kalman filter is shown in Fig. 5.6-1. The algorithm is 
recursive, which means that it processes an estimated state 
vector at time t^ to produce as output an estimated state vector 
for the next sampling time t 2 based on the radar measurements 
at time t 2 * As indicated in the upper left corner of Fig. 5.6-1 
the filter uses the estimated state vector at time t^ to predict 
the states at time t 2 *. This prediction is optimal and takes 
into account the physics of the payload dynamics. 
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STATE VECTOR PREDICTION 


Figure 5.6-1 Block Diagram of Data Processing in a 

Kalman Filter 


Next, the filter uses the predicted state to predict 
what the radar measurements will be at time t 2 * This prediction 
takes into account the radar error model, which represents 
both noise and systematic measurement errors. 


The differences between the actual residual radar 
measurements at time t 2 and the predicted measurements are 
computed. These differences are the one-step-ahead prediction 
errors of the filter and are known as the innovations. Based 
on these innovations, the filter computes an optimal update 
that is added to the predicted state vector for time t 2 . This 
update is optimal; it takes into account the modeled radar 
error sources and the expected accuracy of the state-vector 
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prediction. For simplicity, Fig. 5.6-1 does not portray the 
error covariance calculations that are also performed by the 
filter at each time step. 

5.6.2 Kalman Filter Equations 

The Kalman filter processes the sequence of residual 
radar measurements z^, k= 0, 1, 2, .... which are computed by 
the radar-measurements module described in Section 5.5. The 
filter parameters change with time; at time step k they are 
contained in the following matrices (with n = number of state 
variables and m = number of scalar radar measurements at each 
sampling time, e.g., m = 3 for azimuth, elevation, and range 
data from one radar): 


II 

-0- 

state transition matrix 

[nxn] 

Ok = 

state process noise covariance 

[nxn] 

H k = 

measurement matrix 

[mxn ] 

R k = 

measurement noise covariance 

[mxm] 

The filter uses these four matrices, together with an initial 
estimate of the state vector, x Q , and its error covariance 
matrix, P Q , to compute the following seven matrices for k = 0, 
1,2,...: 

x = 

one-step-ahead estimate of the 
state vector 

[nxl] 

P = 

error covariance of x 

[nxn] 

K 

Kalman gain matrix 

[nxm] 

ii 

innovations vector at 
time k 

[ mx 1 ] 

c k = 

innovations covariance matrix 
at time k 

[mxm] 


5-29 



[nxl ] 


Xj c ( + ) = updated (filtered) estimate of 
the state vector at time k 

P k ( + ) = error covariance of x, ( + ) [nxn] 

at time k 

The filter computes these matrices using the follow- 
ing recursive formulas (in which M [mxn], MM [nxm] , N [nxn] 
are work arrays): 


Initial Conditions 



P 

= 

P o 


[nxn] 

(5. 

6-1) 


a . 

X 

— 

5o 


[nxl ] 

(5. 

6-2) 

For 

k = 

0 

to kmax 






M 

- 

*V P 


[mxn] 

(5. 

6-3) 


C k 

= 

M ' H k + R k 


[mxm] 

(5. 

6-4) 


K 

= 

m T r -l 
M -C k 


[nxm] 

(5. 

6-5) 


N 

= 

P - K-M 


[nxn] 

(5. 

6-6) 


MM 

- 

N-hT 

k 


[nxm] 

(5. 

6-7) 


P k 

( + ) 

= [n - mm-k t ; 

1 + K-R k -K T 

[nxn] 

(5. 

6-8) 


^k 

- 

-k * H k‘- 


[mxl ] 

(5. 

6-9) 


-k 

( + ) 

= x + K-i/ k 


[nxl] 

(5. 

6-10) 


A 

X 

= 

V*k < + > 


[nxl] 

(5. 

6-11) 


p 

- 

v p k <+)- 'l'k + 

<>k 

[nxn] 

(5. 

6-12) 

Next k 



[End of For 

-Next Loop] 



In 

Eqs , 

5. 

6-1 to 5.6-12, 

the following 

calculations 

are 


performed : 
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• The initial state estimate and its error 
covariance are used to initialize x and 
P in Eqs. 5.6-1 and 5.6-2 

• The innovations covariance and Kalman 

gain matrix K are computed in Eq. 5. 6- A 
and Eq . 5 . 6 - 5 

• The error covariance matrix P^( + ) of the 

filtered state estimate is computed 

(using the Josephson-Bierman update) 
(Refs. 5 and 10) in Eqs. 5.6-6 to 5.6-8 

• The innovations vector v, is computed in 

Eq. 5.6-9 K 

• The filtered state estimate x, ( + ) is 
computed in Eq. 5.6-10 

• The one-step-ahead prediction of the 
state vector x and its error covariance 
P are computed in Eqs. 5.6-11 to 5.2-12. 


5.6.3 State-Space Model 


The algorithm represented by Eqs. 5.6-1 to 5.6-12 is 
the optimal filter for estimating the state x^ of the following 
stochastic model for (1) the payload perturbations away from 
the nominal trajectory and (2) the noisy radar tracking data : 


-k+1 = V*k + H k 

[nxl ] 

(5.6-13) 

£ k * « k -x k + v k 

[raxl ] 

(5.6-1A) 

E(w k ) = 0 

[nxl ] 

(5.6-15) 

E[y k J = 0 

[mxl ] 

(5.6-16) 

= Q k '6 k -j 

[ nxn ] 

(5.6-17) 

E!v k -vj) = R k -6 k . j 

[mxm] 

(5.6-18) 
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Elw^vJ) = 0 

[nxm] 

(5.6-19) 

for all k and j 



E l£ k *Wj] = 0 

[nxn] 

(5.6-20) 

for all j > k 



E |X Q ] = 0 

(nxl ] 

(5.6-21) 

ElV^l = P o 

[nxn] 

(5.6-22) 


Equations 5.6-13 to 5.6-22 have the following inter- 
pretation: 


• The state vector x^ is the solution of 

the difference Eq. 5.6-13. The initial 

state vector x is a random variable with 
— o 

zero mean (Eq. 5.6-21) and covariance 
matrix P Q (Eq. 5.6-22). Equation 5.6-13 

is driven by a white noise vector w^. 

• The dimension n of the state vector is 
nine for the filters implemented in this 
study. Six of the states represent the 
dynamical state of the payload relative 
to the nominal trajectory. The remaining 
three states represent radar measurement 
biases or ramps caused by radar calibra- 
tion error. If a bias and a ramp are 
modeled in each of the three measurement 
channels then n=12. 


5.6.4 State Variables 


The state variables (states) are the elements of the 
state vector x^ [nxl]. The first six states are the position 
r^ and velocity v^ of the payload relative to the. nominal tra- 
jectory (expressed in earth-centered Cartesian inertial coordi 
nates). These states form the following 6x1 vector of payload 
states : 
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(payload) 

-k 



[6x1] (5.6-23) 


The remaining states are used to model radar measure- 
ment biases or ramps caused by radar calibration errors. In 
this study three radar error states were used to filter and 
smooth the tracking data provided by WFF from Peruvian radars 
Nos. 8 and 41. All states used to model radar errors are 
elements of the following vector of radar error states (in 
this example there are three error states): 


( radar) 
-k 


e k u) 

e k < 2) 

e k< 3 ) 


[3x1] (5.6-24) 


The complete state vector is organized with the 
payload states listed first: 


(payload) 

-k 

(radar) 

-k 


[9x1] (5.6-25) 


The partitioning of the state vector in Eq. 5.6-25 
induces a partitioning of the parameter matrices 0^, Q^, and 
as indicated in the following equations: 


(payload) 


^k [6x6] 

0 [6x3] 

° 13x6] 

. ( radar) 
^k [3x3] 


[9x9] (5.6-26) 
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(payload) 



[6x6] ! 0 (6x3] 


I n (radar) 
[3x6] • % [3x3] 



(payload) 
k [3x6] 


( radar) 
H k [3x3] 


[9x9] (5.6-27) 


[3x9] (5.6-28) 


The definitions of these partitions and the measurement noise 
covariance matrix are presented in the following sections. 


5.6.5 Transition Matrix 4> 


As indicated in Eq. 5.6-26, the transition matrix is 
block diagonal. This structure occurs because the payload 
states are independent of radar errors. The payload transition 
matrix is partitioned into 3x3 submatrices: 


^ (payload) 
k 


4>n 4> 12 
4>2i $22 


[6x6] (5.6-29) 


The 3x3 submatrices in Eq. 5.6-29 are computed using the algo- 
rithm specified in Section 5.4.3, Eqs . 5.4-50 to 5.4-52 (the 
payload state vector x£P a ^ oac ^ in Eq. 5.6-25 plays the role 
of the state correction 6x^ in Section 5.4.3). 


The transition matrix for the radar error states is a 
3x3 identity matrix when there are three error states: 


^(radar) 


10 0 
0 10 
0 0 1 


[3x3] (5.6-30) 


Equation 5.6-30 is an appropriate transition matrix for modeling 
systematic (bias and ramp) errors in the radar measurements. 

In general, if there are n g (bias and ramp) error states, then 
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/ r^Hpr \ 

<}>k is the n xn identity matrix. (For example, if both 

bias and ramp errors are modeled for azimuth, elevation, and 
range channels, then n = 6.) 

5.6.6 State Noise Covariance Q 

According to the stochastic model used in this study, 
uncertainty in the payload state at time step k is the result 
of uncertainty in the payload state at the initial time k=0. 
Therefore, the state-space model for the payload states is a 
deterministic difference equation with random initial conditions. 
Because there is no white noise driving the payload states, 
the payload Q-matrix in Eq. 5.6-27 is zero: 

q^payload) _ zero matrix [6x6] (5.6-31) 

The radar error states also satisfy a deterministic 
difference equation with random initial conditions. The uncer- 
tainty at the initial time corresponds to uncertainty about 
the radar calibration errors. The radar errors are modeled as 
systematic. Therefore, there is no white noise driving the 
radar error states, and the radar Q-raatrix in Eq. 5.6-27 is 
also zero: 

Q^ rac * ar ^ = zero matrix [3x3] (5.6-32) 

5.6.7 Measurement Matrix H 

The payload measurement matrix H^P a yl° ac *) in £q. 5.6-28 
models linearly the residual tracking measurements that would 
be acquired by an ideal radar. According to the analysis in 
Appendix E, the payload measurement matrix is computed using 
the following equation: 
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H ( payload) _ I" H o 
k I k 

In Eq. 5.6-33, 0^ is the 3x3 zero matrix, and the 3x3 matrix 
is computed using the nominal radar measurements Z^. (The 
algorithm for computing Z k is presented in Section 5.5.) The 
matrix is computed as follows: 

H® = Aj^-B-C k [3x3] (5.6-34) 

In Eq. 5.6-34, the matrices and C k depend on the 
nominal radar measurements (AZ k = azimuth, EL^ = elevation, 
and RA^ = range) as follows: 

Z k = [AZ k EL k RA k ] T [3x1] (5.6-35) 


0 


3 


[3x6] 


(5.6-33) 


Matrix A k is computed using the 

following definitions: 

CL 

cos(EL k ) 



(5.6-36) 

SL 

sin(EL k ) 



(5.6-37) 

CZ 

cos(AZ k ) 



(5.6-38) 

SZ = 

sin(AZ k ) 



(5.6-39) 

RC 

RA k -CL 



(5.6-40) 

RS 

RA k -SL 



(5.6-41) 


~-RC-SZ 

-RS* CZ 

cz-cl" 


ii 

< 

RC-CZ 

-RS* SZ 

SZ-CL 

(5.6-42) 


0 

-RC 

-SL 


In Eq. 5.6-34 

, matrix B 

depends 

on the geodetic latitude LA 

and longitude 

LO of the 

tracking 

; radar as follows: 
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-sin(LA) *sin(LO) cos(LA) 
cos(LO) 0 

-cos(LA) *sin(LO) -sin(LA) 

(5.6-43) 

Matrix in Eq. 5.6-34 depends on the sampling interval fit 
[s] of the radar measurements and the earth's angular velocity 
ft [rad*s ^ ] as follows: 

cos(a k ) sin(a k ) 0 

-sin(a k ) cos(a k ) 0 (5.6-44) 

0 0 1 

a k = ft* 6 1 • k [rad] (5.6-45) 

/ raHp j* ) 

The radar error states x k represent biases or 

ramps in the residual radar measurements. Whether a particular 
error state represents a bias offset or a linear ramp depends 

( T" P H n >* j 

on the radar measurement matrix H k in Eq. 5.6-28. For 

example, in this study three error states were used to process 

test data from WFF . Therefore, x5 rac * ar ^ is a 3x1 vector and 
/ radar) 

H k is a 3x3 matrix in this case. The radar calibration 

errors are modeled as being statistically uncorrelated between 
the azimuth, elevation, and range channels. Therefore, each 
column of H k rac * ar ^ i s filled with zeros except for one entry, 
which is either a 1 or k. The entry is 1 if the error state 
models a bias; it is k if the error state models a ramp. 

As an example, consider the case where all three 
measurement channels have bias errors, and ramps are not being 
modeled. Suppose furthermore, that the first radar error state 




-sin(LA) •cos(LO) 
B = -sin(LO) 

-cos (LA) *cos(L0) 


5-37 


corresponds to the bias in the first measurement channel (azi- 
muth), the second error state corresponds to the second meas- 
urement channel (elevation), the third error state corresponds 
to the third measurement channel (range), then 


„(radar) 

H k 



0 

1 

0 


0 

0 

1 


(5.6-46) 


If on the other hand, each measurement channel error is modeled 
as a ramp, and pure bias offsets are not modeled, then 


„(radar) 

H k 


k 0 0 
0 k 0 
0 0k 


(5.6-47) 


As a third example, suppose that the first measurement 
channel has both bias and ramp errors, the second measurement 
channel has only a bias error, the third measurement channel 
has neither bias nor ramp error, then 


H (radar) 

H k 


1 k 0 
0 0 1 
0 0 0 


(5.6-48) 


According to Eq. 5.6-48, the first error state repre- 
sents a bias in the first measurement channel, the second error 
state represents a ramp in the first measurement channel, and 
the third error state represents a bias in the second measure- 
ment channel. An alternative to Eq. 5.6-48, which assigns the 
three error states differently, is the following: 


1 

0 

0 


k 

0 

0 


H (radar) _ 
H k 


0 

1 

0 


(5.6-49) 



Equations 5.6-48 and 5.6-49 are different ways of representing 
the same radar error model. The two equations simply assign 
different physical meanings to the three error states. 


If the radar data are to be processed with both bias 


and ramp errors being modeled for all three measurement channels, 
then six radar error states are required. In this case j^ ra< * ar ) 
would be a 6x1 vector and H^ rac ^ ar ) would be a 3x6 matrix, which 
could have the following form: 


H ( radar) 
H k 


1 k 0 0 0 0 
0 0 1 k 0 0 
0 0 0 0 1 k 


(5.6-50) 


In Eq. 5.6-50, the first two error states represent the bias 
and ramp errors in the azimuth data, the third and fourth error 
states represent the bias and ramp errors in the elevation 
data, and the fifth and sixth error states represent the bias 
and ramp errors in the range data. 


5.6.8 Measurement Noise Covariance R 


The measurement noise covariance is a 3 X 3 diagonal 
matrix. This models the noise signals in the measurement 
channels as being uncorrelated with each other, which is con- 
sistent with the spectral coherence plots for radar data studied 
under Task 1 of this investigation. The covariance R^ repre- 
sents the level of white noise in each channel of the residual 
radar measurement vector z^. By appropriately defining R^ as 
a function time (k), the Kalman filter will optimally process 
radar data having time-varying noise and gaps caused by missing 
observations. 
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There are three modes in which the Kalman filter is 
used to process residual radar measurements: 


Mode 1 


Mode 3 


The noise covariance R^ is constant (i.e., 


” ^baseline 


for k = 0, 1, 2 k max^ 


Mode 2 The noise variances in are intentionally 

set to very large values, R^ = R^g, for- 

pre-specified values of k that correspond 
to known intervals of missing data or very 
noisy data. 


Individual noise variances in R^ are 

automatically increased (R fc > R b aseline ) 

by the filtering algorithm when an 
individual azimuth, elevation, or range 
measurement is an outlier ( i . e . , the 
measurement is statistically inconsistent 
with the noise and error models for~ 
which the filter is optimized) . The 
outliers are detected automatically using 
the innovations data generated by the 
filter. The amount of increase in R^ is 

a continuous function of the magnitude 
of discrepancy between the observed inno- 
vation values and their theoretical rms 
values . 


The algorithm can be switched between these three 
modes at any values of k during the processing of the radar 
data. The equations for computing the R^ matrix are presented 
in the following. 

Mode 1 for Fixed Noise Model - The baseline value for 
the measurement noise covariance is used for all k: 


baseline 


0 0 


(5.6-51) 
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Ri — R, --I . 

k baseline 


for all k 


(5.6-52) 


In Eq. 5.6-51, the sigmas are the nominal rms values of the 
white noise in the azimuth, elevation, and range data. The 
baseline covariance matrix is diagonal because the noise is 
modeled as uncorrelated between the measurement channels. 


Mode 2 for Data Gaps - Very large noise variances 
(R^^g) are used for pre-specified values of k and pre-specified 


measurement channels 


R 


big 


of. 0 

big az 


0 

0 


r big el 


0 

0 


0 


big ra 


(5.6-53) 


In Eq. 5.6-53, the sigmas are large (e.g. , 1000 times larger 
than the baseline sigmas in Eq. 5.6-51) for each measurement 
channel that has missing or very noisy data. For example, if 
the range data are missing or very noisy, but the angle meas- 
urements are normal, then a reasonable choice is a. 


1000 • a , while a K . „ „ = a and o, . , = a 

ra’ big az az big el el 


big ra 


R^ = ^big ^ Gr P res P ec if ied k 


(5.6-54) 


Mode 3 for Automatic Outlier Processing - Each diagonal 
element of R^ is computed, on the basis of the innovations 
and their variances. The Kalman filter algorithm computes the 
3x1 vector and its covariance matrix as specified in Sec- 
tion 5.6.2. The elements of are arranged in the same order as 
the elements of the residual radar measurement vector z^ 

(1st = azimuth, 2nd = elevation, 3rd = range). In this dis- 
cussion, the following notation is used: 
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(5.6-55) 


-k = * Z 1 z 2 z 3* 


— k = I'l v 2 v 3 ] 


(5.6-56) 


c k = 0 


C 1 0 


(5.6-57) 


baseline 


R k = 


r 1 (bl) 0 0 

0 r 2 (bl) 0 

0 0 r 3 (bl) 


(5.6-58) 


(5.6-59) 


When Mode 3 is in effect, the diagonal elements of R^ 
are computed based on the squares of the innovations and their 
theoretical expected values. For these computations, the vari- 
ance function V is defined for i = 1, 2, and 3: 


V[^ i ,c i ,r.(bl)) = 


r i (bl) + 4-x* [x/A] 
” 1 + |x/A] 7 


(5.6-60) 


where. 


A - 12* c . 

l 


x = v? 
i 


(5.6-61) 


(5.6-62) 


Given: (1) the baseline noise covariance R, , . : 

baseline’ 

(2) the innovations vector and (3) the innovations covari- 
ance matrix C^; the covariance matrix is computed using the 
following algorithm: 
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For i = 1 to 3 


If x/A > 10 then r. = 4*x (avoids possible 

overflow in 

function V) (5.6-63) 

Else r. = V[y. ,c.,r.(bl)]: ' (5.6-64) 

Next i 

This algorithm for computing matrix yields a model 

noise variance r. that equals the baseline value when the squared 

2 1 

innovation is small compared to 12*c.. When the innovations 
start exceeding approximately three standard deviations, then 
the noise variance r^ starts to approach 4*^. The result of 
this is that large innovations (which are improbable under the 
baseline noise model) cause the filter noise model to be changed 
automatically so that noise spikes in the data are filtered 
optimally with respect to the increased noise variance. 

5.6.9 Kalman Filter Outputs 

Figure 5.6-2 depicts the inputs and outputs of the 
Kalman filter module for the case in which there are nine state 
variables: three payload position states, three payload veloc- 

ity states; and three radar error states. Sixty scalars are 
stored per time step: nine filtered state estimates; 45 distinct 

elements of the state error covariance matrix (symmetric 9x9); 
and six scalars representing the innovations and their variances. 

The outputs of the filter module^are stored in a random 
access data store for subsequent ■ processing by the smoothing 
module. The data store is random access, rather than sequential 
access, because the smoother processes the data backwards in 
time. More specifically, the data store should support last-in 
first-out data accesses. 
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A-3770* 



Figure 5.6-2 Kalman Filter Module 

5.7 SMOOTHER MODULE 
5.7.1 Introduction 

The smoother is the last stage of the trajectory 
estimation algorithm. As indicated in Fig. 5.7-1, the inputs 
to the smoother module are taken from the data store in which 
the outputs from the filter module were saved. This data 
store is used in a last-in first-out mode because the~ smoothing 
algorithm processes the data backwards in time. The inputs to 
the smoother at time step k are the processing mode and the 
following time series (in this example there are three radar 
error states ) : 

x^(+) = filtered state-vector estimate 
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[9x1] (5.7-1) 




A 3769 


SMOOTHED ESTIMATE 
OF STATE VECTOR 


ERROR COVARIANCE 
MATRIX 


Figure 5.7-1 Smoother Module 


p k< + > 

= filtered state error covariance 

( 9x9 ) 

(5.7-2) 

^k 

= innovations 

vector 

[3xl] 

(5.7-3) 

c k 

= innovations 

covariance matrix 

[3x3] 

(5.7-4) 


(only 3 diagonal elements are 
used during mode 3 outlier 
processing) 

These inputs are defined in Section 5.6.2, and the processing' 
modes for missing data and outlier processing are discussed in 
Section 5.6.8. 

The outputs of the smoothing algorithm at each time 
step k are the smoothed estimate of. the state vector and its 
error .covariance: 
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x^ = smoothed state-vector estimate [9x1] (5.7-5) 

P, = smoothed state error covariance [9x9] (5.7-6) 

matrix 


In Eq. 5.7-5, the first six elements of the state vector are 
the smoothed position and velocity coordinates of the payload 
relative to the nominal trajectory X^ defined in Section 5.4. 
(All coordinates are expressed in earth-centered Cartesian 
inertial coordinates.) Therefore, the smoothed estimate of 
the payload position and velocity relative to the center of 
the earth is computed by adding the first six elements of the 
smoothed state in Eq. 5.7-5 to the nominal trajectory de- 
fined by Eq. 5.4-56. The mean-square accuracy of this estimate 
is represented by the error covariance matrix in Eq. 5.7-6. 


The remaining states (there would typically be three 
to six of them) in the state vector of Eq. 5.7-5 are radar 
error states. Their covariances are also contained in the 
covariance matrix P^. 


5.7.2 Smoothing Equations 


The algorithm presented in this section is the Frazer- 
Bryson smoother (Ref. 11). It computes trajectory estimates 
that are fully optimal (unbiased and minimum-variance ) with 
respect to the stochastic model for radar errors and payload 
initial state uncertainty. The smoothing algorithm is mathe- 
matically equivalent to batch processing optimally all of the 
available data. 


The algorithm operates on the input data backwards in 
time, working from the final time step k = k back to the 
initial step k = 0. The equations describing the algorithm 
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are presented in the following (A [nxn], A [nxn], b [nxl], b 

— o — ■ —i 

[nxl], B I nxn], C [nxm] , and S -(nxn) are work arrays): 

Initial Conditions 

b = zero vector [nxl] (5.7-7) 


A = 

: zero matrix 

[nxn] 

(5.7-8) 

For k = 

k max to 0 ste PP in g by - 1 



-k 

= x k <+) - P k ( + )‘b 

[nxl] 

(5.7-9) 

P k 

= P k <+) - P k (+)-A-P k (+) 

[nxn] 

(5.7-10) 

C 

= »w 

[ nxm ] 

(5.7-11) 

S 

- c ; H k 

[nxn] 

(5.7-12) 

B 

= (I n - P k< + >- S )->k 

[nxn] 

(5.7-13) 

b 

— o 

= B T -|b - C- ik l 

[nxl ] 

(5.7-14) 

b 


[nxl ] 

(5.7-15) 

A 

o 


[nxn ] 

(5.7-16) 

' A 

iA o 

[nxn ] 

(5.7-17) 

Next k 




H^, and 

Sections 

urement 

In Eqs. 5.7-7 to 5.7-17, the parameter matrices 4> k , 
R k are computed by using the algorithms specified in 
5.6.5, 5.6.7, and 5.6.8. When computing the meas- 
noise covariance R k> the smoother should use the same 


5-47 


processing mode (mode 1, 2, or 3 as defined in Section 5.6.8) 
that is used in the Kalman filter module. For example, if the 
filter changes from mode 1 to mode 2 at k = 55, then the smoother 
module should compute using mode 1 for k < 55 and switch to 

mode 2 processing at k = 55. Using inconsistent modes in the 
filter and smoother modules can produce covariance matrices 
having negative elements along their diagonals. 

The innovation variances are used only for mode 3 
processing. Therefore, when mode 1 or mode 2 processing is 
used at time k, the innovation covariance matrix is not 
used by the smoother and is not a required input quantity. 


5.7.3 Outputs 

The outputs of the smoother module- ax each time step k 
are the estimated state vector and its error covariance. The 
state vector may be written in the following partitioned form, 
in which the payload states are distinguished from the radar 
error states: 


- (payload) 
-k 

-(radar) 

-k 


Inxl] (5.7-18) 


In Eq. 5.7-18, the payload state vector is [6x1], and the radar 

error state vector is [n xl) where n is number of radar error 

e e 

states in the model. (For processing the WFF radar data from 
Peru, n e = 3.) The partition in Eq. 5.7-18 induces the follow- 
ing partitions in the error covariance matrix: 


p(payload) P (p/r) 
k *k 


P (r/p) p ( radar) 

*k *k 


[nxn] (5.7-19) 
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p^P a yl° a d) _ error cov a ri a nce of ^(payload) [6x6] (5.7-20) 

p^radar) =' error covariance of ££ rac * ar ) [n e xn g ] (5.7-21) 

T 

P^ P//r ^ = p£ r//p ^ = error cross-covari a nces (5.7-22) 

The estimated payload state in Eq. 5.7-18 contains 
the estimated position and velocity of the payload relative to 
the nominal trajectory defined in Section 5.4.3. The posi- 
tion coordinates and velocity coordinates are expressed in 
earth-centered Cartesian inertial coordinates, with the posi- 
tion coordinates listed first: 

[6x1] (5.7-23) ^ 


[6x1] (5.7-24) 

To compute the estimated payload position and velocity 
coordinates relative to the center of the earth, r^ and v^, 
the estimated state in Eq. 5.7-23 is added to the nominal state 
in Eq. 5.7-24: 


- 

Ek ’ tk * 5 k 

[3xl] 

(5.7-25) 

. 

. . = 2 k + v k 

.13x1] 

(5.7-26) 


The error covariances of r^ and v^ are given by matrix p£P a yl° a( *) 
in Eq. 5.7-19. 

The estimated systematic radar measurement errors 
(which are modeled by the radar error . states ) are computed 
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using the radar measurement matrix u^ rac * ar ^. This 3*n matrix 

k € 

is defined by Eq. 5.6-28 and the discussion in Section 5.6.7, 
The estimated bias and ramp errors in the radar measurements 
at time step k are computed as follows: 


(bias/ramp) _ „(radar) -(radar) 
-k k *-k 


(3x1] (5.7-27) 


The error covariance matrix of the estimated bias and ramp 
errors is computed using the following formula: 


cov 



T 

u (radar) -(radar) u (radar) 

H k >P k * H k 

(3x3] (5.7-28) 


These equations are applied to the analysis of WFF radar data 
from Peru in Section 5.8. 


5.8 VERIFICATION OF THE TRAJECTORY ESTIMATION ALGORITHM 

To verify the performance of the trajectory estimation 
algorithm, WFF provided tracking data from Radars No. 8 and 
No. 41 in Peru. The results of processing these data to esti- 
mate a Nike-Orion trajectory and a Terrier-Malemute trajectory 
are discussed in Sections 5.8.1 and 5.8.2. 


5.8.1 



The data used in this verification of the trajectory 
estimation algorithm were obtained from Radars No. 8 and No. 41, 
which were simultaneously tracking a Nike-Orion (31.027) trajec- 
tory. The results of processing data from Radar No. 8 are 
discussed first. Then the results obtained using data from 
Radar No. 41 are presented. Section 5.8.1 concludes with a 
comparison of the two independent estimates of the trajectory. 
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Radar No. 8 Data - Figure 5.8-1 depicts the residual 

tracking data from Radar No. 8 used for this test. These plots 

show the departures of the actual measurements from the nominal 

measurements. (The nominal measurements would have been obtained 

if an ideal radar had tracked the nominal trajectory.) The time 

is measured starting with the nominal apogee, and the initial 

conditions for the nominal trajectory are consistent with the 

position and velocity data provided in the WFF documentation. 

The parameters for the atmospheric drag force model are as 

follows: drag coefficient = 1.5; payload mass = 286 kg; and 

2 

payload cross-sectional area = 0.15 m . 




Figure 5.8-1 Residual Tracking Data, Nike-Orion (31.027) 

Trajectory, Radar No. 8 
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For these validation tests, the radar data (originally 
sampled at 1-s intervals) were intentionally undersampled with 
a sampling time of 5 s between consecutive measurements. This 
provided. a realistic test of the trajectory estimation algorithm 
while using a reduced number of measurements. 

The error model for Radar No. 8 had three state vari- 
ables: one of the states modeled a bias in the azimuth channel; 
while the other two states modeled a bias and a ramp in the 
range channel. The rms a priori uncertainties of these errors 
were 0.5 deg for the azimuth bias, 1000 m for the range bias, 
and 6 m/s for the range ramp. In addition to these calibration 
error uncertainties, each measurement channel was modeled as 
having random white noise with the following rms values: 

0.04 deg in azimuth, 0.09 deg in elevation, and 3.7 m in range. 
These noise values are the rms of residual tracking data com- 
puted from raw measurements by subtracting a least-squares 
linear trend from each channel of data. For estimating the 
noise levels, typically 50 to 100 samples of radar data sampled 
at 1-s intervals were used. 

The a priori rms uncertainties of the initial payload 
position and velocity were 1000 m for each of the three position 
coordinates and 10 m/s for each of the velocity coordinates. 
These values were selected to be large but plausible so that 
the filter/smoother would rely primarily on the tracking data, 
rather than on the initial conditions, for estimating the 
trajectory. 

The behavior of the Kalman filter can be monitored by 
observing the innovations data, which are the filter’s one-step- 
ahead predictions of the residual radar data. Figure 5.8-2 
depicts the innovations for the Nike-Orion data. In these 
plots the innovations have been normalized (divided by their 
theoretical standard deviations) so that outliers can be more 
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Figure 5.8-2 Kalman Filter Innovations, Nike-Orion (31.027) 

Trajectory, Radar No. 8 


easily identified. For example, in the plot of the range 
innovations, there is a 4-sigma outlier (a measurement that 
lies 4 standard deviations away from its expected value of 
zero) . The theoretical standard deviations of the innovations 
are the square roots of the diagonal elements of the innovations 
covariance matrix C^. The covariance is computed automat- 
ically by the Kalman filter using Eq. 5.6-4. When the stochastic 
error model of the Kalman filter is consistent with the tracking 
data, the innovations are samples of unit-variance zero-mean 
white noise. Statistically significant departures from this 
behavior indicate that the data may contain errors that are not 
modeled by the filter. 
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Figure 5.8-3 shows the outputs of the smoother module. 
In the first column of Fig. 5.8-3, the position of the payload 
relative to the nominal trajectory is plotted as a function of 
time. The solid lines are the smoothed estimates, and the 
dotted lines are the 1-sigma uncertainties (one-standard-devia- 
tion error bounds) of these estimates. The estimated payload 
velocity relative to the nominal trajectory is plotted on the 
right side of Fig. 5.8-3. Both position and velocity are 
expressed in earth-centered Cartesian inertial coordinates. 


RESIDUAL POSITION 
ESTIMATES 



TIME lied 



RESIDUAL VELOCITY 
ESTIMATES 



TIME ls«cl 


Figure 5.8-3 Smoothed Trajectory Estimates, Nike-Orion 

(31.027) Trajectory, Radar No. 8 
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The following comments apply to Fig. 5.8-3: 


• The rms position accuracy ranges from 
110 m to 540 m, depending on the position 
coordinate and the time 

• The rms velocity accuracy ranges from 
1.2 m/s to 3.4 m/s with very little 
dependence on the time 

• The rms accuracy estimates are computed 
by the algorithm, based on the error 
covariances generated by the Kalman 
filter/smoother 

• The algorithm was processing the data in 
mode 3 for automatic outlier detection 
and handling. Therefore, the outlier in 
the range innovations was automatically 
detected and appropriately processed. 

The algorithm also estimated the systematic bias and 
ramp errors in the radar measurements. The smoothed estimates 
(o = theoretical standard deviation of the estimation error) 
are as follows: 


0.05 deg (o = 0.46 deg) 
1.8 m/s (ct = 3.5 m/s) 

200 ra (ct = 470 m) 


Azimuth Bias = 
Range Ramp = 
Range Bias = 


These results indicate that the bias and ramp estimates are 
imprecise (because the standard deviations are larger than the 
estimates of the biases and the ramp). The reason for this is 
that data from a single radar provides insufficient information 
for the precise estimation of radar errors. Better precision 
can be obtained by simultaneously processing the tracking data 
from two or more radars tracking a single payload. The Kalman 
filter/smoother can be extended to handle data from multiple 
radars, but such an extension was not part of this investigation. 
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Radar No. 41 Data - The residual data from Radar 
No. 41, tracking the Nike-Orion 31.027 trajectory are depicted 
in Fig. 5.8-4. An extremely large outlier occurs at 15 s in 
the range data. This outlier is caused by 3 data processing 
error and would normally be edited manually. However, as a 
demonstration of the algorithm's automatic outlier processing 
(mode 3 processing), the outlier is intentionally left in the 
data. As with Radar No. 8 , the data in Fig. 5.8-4 are sampled 
at 5-s intervals. 





Figure 5.8-4 Residual Tracking Data, Nike-Orion 

(31.027) Trajectory, Radar No. 41 
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Three error states were used to model systematic bias 
errors in the three radar channels. The a priori rms uncer- 
tainties of these errors were 0.5 deg for azimuth and elevation 
and 1000 m for range. In addition the random errors .in each 
channel were modeled as white noise with the following rms 
values: 0.04 deg in azimuth; 0.03 deg in elevation; and 45 m 

in range . 


The innovations from the Kalman filter are depicted 
in Fig. 5.8-5. The data are normalized by their theoretical 
standard deviations. Except for the isolated outlier in the 
range data, the innovation values lie in the expected range. 






Figure 5.8-5 Kalman Filter Innovations, Nike-Orion 

(31.027) Trajectory, Radar No. 41 
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The smoothed position and velocity estimates, and 
their 1-sigma error bounds, are presented in Fig. 5.8-6. The 
following comments apply to these results: 


• The rms position accuracy is in the range 
from 370 m to 580 m 

• The rms velocity accuracy is in the range 
from 1.6 m/s to 2.2 m/s 

• The algorithm was operating in mode 3 
and therefore automatically ignored the 
extreme outlier in the range data. 
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Figure 5.8-6 Smoothed Trajectory Estimates, Nike-Orion 

(31.027) Trajectory, Radar No. 41 
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The smoother also provided the following estimates of 
the radar biases (o = theoretical standard deviation of the 
estimation error): 

Azimuth Bias = 0.05 deg (a = 0.46 deg) 

Elevation Bias = 0.04 deg (a = 0.18 deg) 

Range Bias = 220 m (a = 420 m) 

The conclusion to be drawn from these results is that the data 
provided by Radar No. 41 alone are not sufficient for precise 
estimates of the radar biases. 

Comparison of Trajectory Estimates - The data from 
Radars No. 8 and No. 41 were processed separately to provide 
two independent estimates of the Nike-Orion trajectory. In 
this section, the two trajectory estimates are compared with 
each other. The comparison shows that the two trajectory esti- 
mates are statistically consistent with each other . 

Table 5.8-1 lists the mean differences between the 
two trajectory estimates (expressed as individual position and 
velocity coordinates averaged over the span of the data proc- 
essed). If the trajectory estimates from the two radars are 
consistent, then the mean differences between them should be 
less than, say, two standard deviations of the difference. 

The theoretical standard deviation (a) of each mean difference 
is also listed in Table 5.8-1. In the last column of the table, 
the normalized mean differences are listed. These normalized 
quantities express each difference as a multiple of its theo- 
retical standard deviation. For consistency at a 2 -ct level, 
the numbers in the last column of Table 5.8-1 should be less 
than 2. An examination of the table shows that this criterion 
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TABLE 5.8-1 


STATISTICS OF DIFFERENCES BETWEEN 
TWO NIKE-ORION <31.027) TRAJECTORY ESTIMATES 


ESTIMATED 

QUANTITY 

MEAN 

DIFFERENCE 

THEORETICAL 
STD. DEV. (a) 

NORMALIZED 
MEAN DIFFERENCE 



1st Position 
Coordinate 

300 m 

410 m 

0.7 a 

2nd Position 
Coordinate 

730 m 

550 m 

1.3 a 

3rd Position 
Coordinate 

620 m 

620 m 

1.0 o 

1st Velocity 
Coordinate 

1 . 6 m/s • 

2 . 1 m/s 

0.8 a 

2nd Velocity 
Coordinate 

3 . 7 m/s 

3.6 m/s 

1.0 a 

3rd Velocity 
Coordinate 

1 . 1 m/s 

3 . 1 m/s 

0.4 a 


for consistency is easily met; the largest normalized differ- 
ence between the two trajectory estimates is only 1.3 a. The 
conclusion to be drawn from this comparison is that the trajec- 
tories estimated from the two radar data sets are statistically 
consistent . 

5.8.2 Terrier-Malemute Trajectory 

The final set of test data used for validating the 
trajectory estimation algorithm was obtained from Radar No. 8 
tracking a Terrier-Malemute (29.019) trajectory. An initial anal 
ysis of the raw tracking data disclosed a significant growing 
oscillation in the range rate having a frequency of 0.5 Hz and 
a peak magnitude of 300 m/s. The cause of this oscillation is 
unknown . 
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Therefore, the range data were smoothed by computing the 2-s 
running mean of the data to suppress the 0.5-Hz error signal. 
The radar data were then resampled at 10-s intervals for proc- 
essing with the trajectory estimation algorithm. In an oper- 
ational setting this resampling would not necessarily be 
recommended. But for validating the trajectory estimation 
algorithm, the resampling is advantageous because it permits a 
realistic test of the algorithm with a smaller data processing 
effort . 


WFF data on the nominal velocity at time k = 0 (nominal 
apogee) were not available for this analysis. Therefore, the 
inital conditions for the nominal trajectory were determined 
by iteration as explained in the following paragraphs. 

The raw radar data (sampled at 1-s intervals) were 
first transformed to earth-centered Cartesian inertial coordi- 
nates. Second, linear trends were fitted to the first 20 s of 
the data. Finally, the linear trends were used to estimate 
the payload velocity coordinates at time step k = 0. The nominal 
trajectory was computed using this estimate of the initial 
payload velocity together with the nominal payload position 
determined from WFF documentation. The nominal measurements 
and residual measurements were computed using the nominal tra- 
jectory, and the residual measurements were filtered and 
smoothed. The a priori state error covariance was the same 
one used for processing data from Radar No. 41 on the Nike- 
Orion trajectory, and the radar error model had three bias 
states (one for each measurement channel) and additive white 
noise. The initial rms uncertainties of the bias states were 
0.5 deg in angle and 1000 m in range. The rms of the white 
noise for each channel was determined using the method recom- 
mended in Section 4.3.2 (based on the level at high frequencies 
of the estimated power spectrum for each data channel). The 
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filter and smoother were operated in mode 1 (no outlier detec- 
tion) so that the measurement noise model had a fixed covariance 
matrix . It is important to use mode 1 processing (and mode 2 
processing if there are gaps caused by missing data) because 
large innovations may be produced during the first iterations, 
when the nominal trajectory is inaccurate. In mode 3 processing, 
these large innovations would be interpreted as outliers and 
the speed of convergence of the iterative process would be 
reduced . 


The smoothed estimate of the initial payload vector 
was then used to update the original initial conditions for 
computing the nominal trajectory. Using the updated initial 
position and velocity, a new nominal trajectory was computed, 
and the whole process was repeated using the same a priori ■ 
state error covariance and measurement noise covariance. The 
smoothed estimates of the state vector were significantly 
smaller than in the previous iteration, which indicated that 
the second nominal trajectory was closer to the true trajectory 
than the first one. 

The resulting smoothed estimate of the payload's ini- 
tial state was again used to update the initial position and 
velocity to compute a third nominal trajectory, and the data 
analysis was repeated a third time. However, this time the 
covariance matrix of the measurement noise model was adjusted 
so that the innovations would have the expected 1-2 sigma range 
of values. The resulting rms measurement noise model was 0.3 deg 
in azimuth, 0.02 deg in elevation, and 15 m in range. The 
results of this last iteration are discussed in the following. 

The residual tracking data from Radar No. 8 are de- 
picted in Fig. 5.8-7, with a sampling interval of 10 s between 
measurements. To model the offsets in the data, three radar 


5-62 


2.0 


A-3933 


RESIDUAL AZIMUTH 



0.5 1 * 1 — 1 

0 100 200 300 

TIME (sec) 



i i_ i J 

0 100 200 300 


TIME (see) 

Figure 5.8-7 Residual Tracking Data, Terrier-Malemute 

(29.019) Trajectory, Radar No. 8 


error states were used to model possible biases in each of the 
three channels. The innovations produced by the Kalman filter 
operating in mode 1 are shown in Fig. 5.8-8. Three-sigma out- 
liers occur in the elevation and range data near the end of 
the data set. Since the filter was operating in mode 1, it 
used a measurement noise model with a fixed covariance matrix. 
Therefore, the outliers were processed suboptimally in this 
example. Nevertheless, the smoothed trajectory estimates de- 
picted in Fig. 5.8-9 are free of jumps and are qualitatively 
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Figure 5.8-8 Kalman Filter Innovations, Terrier-Malemute 

(29.019) Trajectory, Radar No. 8 


similar to the results of mode-3 processing for the Nike-Orion 
trajectory discussed in Section 5.8.1. Smooth trajectory 
estimates are produced in this case because the filter/smoother 
produces a trajectory estimate that is consistent with the 
physics of a massive payload in free fall. 


An examination of Fig. 5.8-9 leads to the following 
conclusions: 
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Figure 5.8-9 


Smoothed Trajectory Estimates, Terrier-Malemute 
(29.019) Trajectory, Radar No. 8 


• The rms position accuracy ranges from 
280 ra to 570 m, depending on the posi- 
tion coordinate and the time 


• The rms velocity accuracy ranges from 
0.6 m/s to 1.7 m/s with very little 
dependence on the time. 


The smoothed estimates of the radar biases are listed 
in the following (o = theoretical standard deviation of the 
estimation error): 
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Azimuth Bias = 0.24 deg (o = 0.40 deg) 

Elevation Bias = -0.24 deg (a = 0.03 deg) 

Range Bias = 1040 m (o = 200 m) 

These results indicate that the elevation and range biases 
have been estimated with reasonable precision because the 
one-sigma error bounds are much smaller than the magnitudes of 
the estimated biases. 

5.8.3 Conclusions 

The validation tests described in Section 5.8 verify 
that the trajectory estimation algorithm yields consistent 
results when it is used to estimate the trajectory of a single 
payload that was tracked simultaneously by two radars. The 
validation tests also verify that the algorithm automatically 
handles isolated data outliers and provides estimates of the 
rms accuracies of the estimated position and velocity coordi- 
nates of the payload. The test results include a demonstration 
of using the algorithm iteratively to estimate the initial 
velocity of a payload when only nominal position data were 
available. 


5 . 9 SUMMARY 

An algorithm has been developed for processing radar 
tracking data to estimate payload trajectories above the atmos- 
phere (altitude > 50 kra). The algorithm is based on Kalman 
filtering and smoothing techniques and is optimal with respect 
to models' for gravitation, atmospheric drag, radar measurement 
errors, and errors in the assumed initial position and velocity 
of the payload. 
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The algorithm has three operating modes: 


• The first mode employs a radar noise 
model having a fixed noise covariance 
matrix. This mode is used for itera- 
tively refining an initial poor estimate 
of the payload initial position or 
velocity . 

• The second mode employs a variable radar 
noise model having a very large noise 
covariance for specified measurements in 
the tracking data. This mode is used 
for optimally processing data having 
gaps caused by missing data, or very noisy 
data. 

« The third mode employs a variable radar 
noise model which employs a noise covari- 
ance that is automatically increased in 
a continuous fashion when the filter 
innovations indicate data outliers. This 
mode is used for optimally processing 
data having isolated data errors, as 
opposed to intervals of many missing or 
noisy data, which are better processed 
using mode 2. 


The algorithm is optimal in the sense that it computes 
unbiased minimum-variance estimates based on the following 
error models: 


• The initial payload position and velocity 
coordinates of the payload are random 
variables with a specified error covari- 
ance matrix 

• The radar noise is modeled as additive 
white noise with a specified baseline 
covariance matrix. In mode 2 or mode 3 
processing, the noise model covariance 
is increased from the baseline values as 
appropriate 
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• The systematic radar errors are modeled 
as biases and ramps. The initial uncer- 
tainty about the magnitudes of these 
systematic errors is represented by an 
error covariance matrix. 

In addition to the error models, the algorithm also 
uses an accurate model of the dynamics of the payload as it 
falls in the earth's gravitational field subject to nominal 
atmospheric drag forces. The Kalman filter and smoother algo- 
rithms are linear time-varying data processors that process 
residual radar measurements (residual data = raw data - nominal 
data, where the nominal data correspond to an ideal radar 
tracking the nominal payload trajectory based on nominal initial 
conditions, the normal gravitational field of an ellipsoidal 
earth model, and nominal atmospheric drag forces). 

Because the algorithm is optimized with respect to 
the physics of the falling payload and the statistics of the 
important error sources, the estimated trajectories are smooth 
even when the tracking data are noisy. Moreover, the algorithm 
automatically computes the error covariance of the estimated 
payload state and radar error states at each time step. 

Using radar data provided by WFF from Radars No. 8 and 
No. 41 in Peru, the performance of the algorithm has been suc- 
cessfully demonstrated. In particular, the algorithm provided 
two mutually consistent estimates of a Nike-Orion trajectory 
that had been simultaneously tracked by the two radars. 
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6 . 


SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS 
FOR FUTURE STUDY 


6 . 1 SUMMARY 


The principal accomplishments of this study are sum- 
marized in the following: 


• The accuracy of the current WFF data 
smoothing technique was analyzed for a 
variety of radars and payloads, using 
tracking data provided by WFF for this 
study 

• Alternative data noise reduction techniques 
were assessed and recommendations were 
made for improving radar data processing 

at WFF 

• A data-adaptive algorithm, based on Kalman 
filtering and smoothing techniques, was 
developed for estimating payload trajec- 
tories above the atmosphere from noisy 
time-varying radar data 

• The new trajectory estimation algorithm 
was tested and verified using radar 
tracking data provided by WFF. 


6.2 CONCLUSIONS 

The principal conclusions of this. study are summarized 
in the following. More detailed discussions are provided in 
Sections 2. A, 3.7, 4. A, and 5.9. 
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Rms noise levels in smoothed radar data vary with payload . 


• The estimated rms noise levels of posi- 
tional data products (for Zuni , Super 
Loki Optical , and Super Loki Sphere 
trajectories tracked by radars Nos. 3 
and 5) produced by the current WFF 
smoothing filters have the following 
range of values for the data sets 
analyzed in this study: 

Latitude and Longitude: 3 pdeg to 73 pdeg 

Height: 1.3 ft to 6.9 ft 


Current WFF smoothing techniques can be improved by using 
the following recommended techniques : 


• Subtracting orthogonal polynomials from 
the radar data before smoothing, to reduce 
possible distortion of the nominal tra- 
jectory by the smoothing filter 

• Estimating the high-frequency rms noise 
levels in smoothed tracking data using 
an autoregressive modeling algorithm, to 
monitor data quality and provide quanti- 
tative error estimates with the existing 
smoothing filters 

• Replacing the current midpoint smoothing 
filters with an alternative Kalman filter/ 
smoother, to process radar data optimally 
when the data contain significant time- 
varying noise, measurement gaps caused 

by missing data, or large uncertainties 
on radar calibration errors. 


The Kalman filter/smoother algorithm yielded reasonable 
trajectory estimates and error analyses when it was tested 
using radar tracking data provided by WFF. 


• The algorithm is data-adaptive , automati- 
cally handles data outliers, and can 
optimally smooth data sets having measure- 
ment gaps caused by missing data. 
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6.3 


RECOMMENDATIONS FOR FUTURE STUDY 


Based on the results of this study, the following 
recommendations are made for future study: 


The Kalman f ilter/smoother algorithm can be extended to 
increase its flexibility : 


• To process data from two or more radars 
simultaneously 

• To model data errors that are more complex 
than data gaps, additive white noise, 
linear trends, and isolated outliers 

• To estimate automatically the initial 
payload position or velocity from the 
radar tracking data alone 

• To calibrate radars based on future 
Geosat satellite tracking techniques. 


A smart preprocessor for radar tracking data can be 
developed: 


• To select and setup appropriate error 
models automatically for the Kalman 
filter/smoother algorithm 

• The preprocessor could employ artificial 
intelligence techniques together with 
automated stochastic state-space modeling 
algorithms. 
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APPENDIX A 

ORTHOGONAL POLYNOMIALS 


The purpose of this appendix is to describe the rec- 
ommended way of using least squares to fit orthogonal poly- 
nomials to time-series data. In this study, the orthogonal 
polynomials are used to model nominal payload trajectory sig- 
nals in radar tracking data. 

The polynomials used in this study are optimal in the 
following sense (Ref. 12): 



This is a conventional least-squares problem using ordinary 
polynomials. The equations for computing the optimal coef- 
ficients become ill conditioned as N and m increase. This 
leads to serious numerical difficulties in practice. To reduce 
these numerical problems, the problem is reformulated using 
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polynomials f n (x) of degree n that are orthogonal on the dis-r 
crete interval k = 0, 1, N-l and normalized so that their 

values are in the range from -1 to +1. A linear combination 
of these polynomials is used in place of the formula for p n (x) 
given in Eq. A-2. They are called Chebyshev polynomials on a 
discrete domain (Ref. 12) and are computed recursively as 
follows. 

f°(x) = 1 (A-4) 

f 1 (x) = 1 - <2x)/N (A-5) 

and for n = 1, 2, .... N-l 

f n+1 (x) = [ (2n+l )(N-2x)f n (x) - n(N+n+l ) f 11 " 1 (x ) ] (n+1 )“ 1 (N-n ) " 1 

(A- 6 ) 

The optimal polynomial p ra (x) of degree m is expressed in terms 
of the orthogonal polynomials as follows: 

m 

P m (x) = 2 a n fn(x) <A-7) 

n=o 

In Eq. A-7, the coefficients a n are computed as follows: 

a = ^ , n = 0, 1, .... m (A-8) 

n (f n >f n ) 

where the numerator is the inner product 

N-l 

(y,f n ) = 2 y k fn(k) <A-9) 

k=o 

A-2 



and the denominator of Eq. A-8 is given by the following ex- 
pression involving factorials 


( f n f n ) - (N+n+1 ) ! (N-n) ! 

<2n+l)(N!) 2 


(A-10) 
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APPENDIX B 

AUTOGRESSIVE MODELING 

o 

This appendix describes an effective method for using 
autoregressive (AR) modeling to estimate the power spectra of 
time series. The AR models are used in this study to estimate 
the power spectra of noise- like errors in radar tracking data. 

An autoregressive model of order p for time series 
y^, k = 0, 1, ..., N-l, is the difference equation 

P 

y k = ^2 C j' y k-j + w k ’ k = p, p+1, ..., N-l (B-1) 
j=l 

Equation B-1 is driven by the residual noise w^. AR models 
are developed from time series by choosing the coefficients 
c j , j = 1 to p, so that the sample mean square (VAR) of the 
residuals is minimized: 

N-l 

VAR = STp E w k (B - 2> 

k-p 

In Eq. B-2, the limits of the summation are chosen to avoid 
running off the ends of the time series. Choosing the AR coef- 
ficients to minimize VAR is known as the covariance method of 
AR modeling. 

If the AR model is appropriate for the process gener- 
ating the data y^, then the residuals w^ are a sample of approx- 
imately white noise. It follows that the power spectral density 
(power spectrum) of the discrete- time process generating the 
data y^ can be estimated as follows: 

B-1 



S C (F) = 


(B-3 ) 


c k' e 


i2nFk 


where F [cycles/sampling interval] is the normalized frequency 
variable. F = 0.5 is the folding frequency. The variance of 
the random process having the power spectrum S Q (F) is given by 
the area under the spectrum (including negative frequencies): 


variance 


= / S o (F)-dF 

* 1/2 


(B-4) 


A natural estimate for the power spectrum of the under- 
lying continuous- time process (of which the data y^ are 
uniformly-spaced sample values, i.e., y^ = y(k*6t) with 
6t = sampling interval) is 


S o (f/f s> 
S(f) = ° . 5 


where 


f = l/(6t) = sampling frequency [Hz] 
f = spectrum frequency [Hz] 

S(f) = power spectrum [variance/Hz] 

For each time series y^, k = 0, 1, ..., N-l, the best 
choise of the order p is estimated by computing the Akaike 
information criterion (Refs. 2-4) for each model in a family 
of AR models. 
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Each model in the family corresponds to a dif- 
ferent model order p = 0, 1, . .., N/20. For each of these 
models the AIC is computed: 

AIC = N-log e (VAR) + 2-p (B-6) 

That model for which the AIC is smallest is chosen as the best 
AR model in the family for the purpose of modeling the under- 
lying process that generated the observed data y^. An algorithm 
(ACOVAR) for efficiently computing the family of AR models and 
selecting the model order using the AIC is specified in Ref. 13. 
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APPENDIX C 

STOCHASTIC STATE- SPACE MODELING 


C . 1 INTRODUCTION 

This appendix describes a canonical-variates (CV) 
technique for the stochastic modeling of vector time series. 
(This is a modification and extension of ideas originally pre- 
sented in Ref. 14). The CV method differs significantly from 
the one-step linear prediction techniques that are now commonly 
used to develop autoregressive (AR all-pole) and autoregressive 
moving- average (ARMA pole-zero) models from empirical data. 
The technique has been used successfully with a variety of 
geophysical data sets for spectrum estimation, reduced-order 
modeling, and optimal filtering. The CV algorithm is recursive 
in the data. Therefore the data may consist of several short 
time series instead of one contiguous sequence. 

The modeling technique was motivated by Akaike's orig- 
inal work (summarized and extended in Ref. 3), which describes 
an ARMA modeling technique based on a canonical-variates analy- 
sis (Ref. 15) in the time domain. W.E. Larimore modified and 
extended Akaike's approach to state-space modeling (Ref. 16). 
In this appendix the CV technique is interpreted as a form of 
multi-step linear prediction . This leads to a simplified deri- 
vation of the algorithm and shows that the technique does not 
reduce to conventional one-step rediction-error modeling. 

An important advantage of the CV approach is that a 
family of optimal state-space models is generated by solving a 
finite number of linear equations for the state-space model 
parameters. This fact distinguishes the CV technique from 
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other approaches, such as Gaussian maximum likelihood and con- 
ventional one-step prediction, which all require that nonlinear 
equations be solved for state- space modeling. These nonlinear 
equations lead to iterative calculations that may or may not 
converge, depending on the particular data being analyzed. In 
contrast, the only iteration in the CV approach occurs while 
computing Singular Value Decompositions (SVDs). The SVD is 
well understood and can be implemented with dependable algo- 
rithms (Ref. 17). 

C . 2 SUMMARY OF RESULTS 

The CV approach to state-space modeling consists of 
three steps. The first step is to solve a family of least- 
squares multi-step linear prediction problems with different 
rank constraints on the predictors . This is a canonical - 
variates analysis of the joint behavior of the local past and 
local future of the time series, in which several future vec- 
tors of the time series are being simultaneously predicted 
using several past vectors. The entire family of predictors 
is computed at one time by using the SVD. The output of this 
analysis is the definition of the canonical states in terms of 
the observed data. The state vector is expressed as a linear 
combination of the local past, and the individual states are 
canonical variates. 

The second step is to solve a family of linear least- 
squares problems that use the state vectors from Step One as 
given quantitites . This yields the parameter matrices of a 
family of state-space models for the time series. These models 
differ from each other in their complexities, ranging from the 
simplest white-noise model containing no states, to a model of 
maximum complexity containing the largest number of states 
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permitted by the CV analysis. For the number of states it 
contains, each of these models is optimal in a least-squares 
sense. 


The third step is to select one of the state-space 
models from the family generated during Step Two. For th< 


For the 


purpose of modeling the underlying process that generated the 
empirical data, the Akaike information criterion (Refs. 2-4) 


is used to select that model which is best supported by the 
data. Alternatively, a model with a reduced number of states 


is selected for reduced-order modeling. A formula is given in 
Section C.3 for computing the amount of mutual information be- 
tween the future and past that is lost by reducing the number 
of states when modeling Gaussian processes. 


C.3 ANALYSIS 


C.3.1 Step One: Canonical Variates 


A time series of empirical data is denoted by the 
sequence of m-vectors y k for k = 1 to n ' . An ergodic state- 
space model (in innovations form) for the data process is 
represented by the following equations: 


k+1 " ^k + G -k 

(C-l) 

Z k = Hx k + y k 

(C-2 ) 

R = 

(C-3) 


The nxl state vector is 
vector at time step k. 


x^, and y k is the observed m*l output 
The nxn state- transition matrix is <t> , 


while G is the nxm noise-gain matrix, and H is the mxn output 
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matrix. The mxm covariance matrix of the zero-mean white noise 

(steady-state) innovations process ^ is R. The noise vector j/j 

is uncorrelated with v, for all j j* k; is also uncorrelated 

k j 

with x^ for all j ^ k. 


The primary objective of state-space modeling is to 
use observations of y^ for k = 1, 2, ...» n' to estimate val- 
ues for the model order, n, and the parameter matrices 4> , G, 

H, and R. The CV technique for doing this is based on a local 
past and local future for v^. The local past z^(p)» of length 
p, is the pmxl vector containing the p most recent predecessors 
of y k : 


z U <P> = 


*k-l 

*k-2 


Uk-p J 


k = p+1 , n'+l ( C - A) 


The local future z^(f), of length f, is the 
tor containing y^ and the next f-1 data vectors: 


fmxl yec- 


4(f) 


*k 

-k+1 


L z k+f-l J 


k = 1, . . . , n’-f+l ( C - 5 ) 


In many applications the lengths of the local future and past 
are equal (f=p). 


The state vector x^ contains n numbers: the states 

at time k. These states contain all the information available 
from the entire past z^(°° ) about the entire future z^O.)* In 
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other words, x k is a sufficient statistic of z k (<») for opti- 
mally predicting z£(°°). 

The empirical data are the finite sequence y k for 
k = 1 to n ' . To estimate the corresponding sequence of state 
vectors, x k is expressed as a linear transformation L(n) of 
the local past z^(p): 

x k = L(n) z k ( p ) , k = p+1, ..., n' (C-6) 

In Eq. C-6, the n*pm matrix L(n) is to be defined so that three 
conditions are satisified: 


• Given n, the vector x^ is optimal for 

predicting the local future z k (f). That 

is, for some fm*n matrix M(n), the fol- 
lowing estimates of the local future are 
optimal in a weighted least-squares sense: 

z k (f) = M(n) x k , k = p+1, .... n’-f+l 

(C-7 ) 

• The state vector x k is standardized (for 

convenience) so that the n states are 
uncorrelated with each other, and each 
state has zero mean and unit sample vari- 
ance. This means that the sample covari- 
ance matrix of x k (defined in Eq. C-8 

and denoted (x,x)) is the n*n identity 

matrix as indicated in Eq. C-9: 


n ' - f +1 

(x,x) - x R x k • [n'-f-p+l]' 1 (C-8) 

k=p+l 

(x,x) = I n (C-9) 
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• The states in x^ are arranged in the 

order of their importance for predicting 
the future (as measured by a weighted 
least-squares error criterion). The 
first state is most important, the second 
state is second most important, etc. 


The optimal L(n) and M(n) matrices are determined by 
solving the following linear prediction problem . Fix the posi- 
tive integers n, f, and p, with n < min[fm,pm]. Find the n*pm 
matrix L(n) and the fmxn matrix M(n) , such that the following 
estimate of the local future 

z^(p) = M(n) L(n ) zj^(p), k = p+1 to n'-f+l 

(C- 10) 

has the error vector 

e k (n) = z k (f) - z k (f) (C-ll ) 

with the smallest weighted sum-square error J(n): 
n'-f+l 

J(n) " ]£ -k (n) W < f >” 1 £ k <n) (C- 12) 

k=p+l 

The fmxfm weighting matrix W(f) is the sample covariance matrix 
of the local future: 

n’-f+l 

W(f) $ ^ z^(f).[z k ( f )] T -[n'-p-f+l]- 1 (C-13) 

k=p+l 

This choice of the weighting matrix leads to a definition of 
the state vector that maximizes the mutual information between 
the local future and the local past. Moreover, W(f) makes the 
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canonical analysis independent of the units in which the ob- 
served data are expressed. (In Eq. C-13), if W(f) is sin- 
gular, then f is larger than it needs to be, and its value 
should be reduced.) 

Fragmented data sets, consisting of several time ser- 
ies, are handled in Eqs. C- 12 and C-13 by summing over all 
contiguous data segments, taking care to avoid running off the 
ends of the segments, and then dividing by the total number of 
terms in each sum. 


The optimal choices for L(n) and M(n) are found by 
defining the standardized future vectors (with their depend- 
encies on p and f suppressed): 

- <z + ,zV T/2 • z£, k = p+1, ..., n'-f+l 

(C-1A) 

In Eq. C-14, the notation (z ,z ) ' represents the trans- 

posed inverse of any matrix square root of the sample covari- 
ance of the local future. The covariance matrix and its 
square-root matrix are defined by the following equations: 

n 1 - f+1 

U + ,z + ) - V z£ [z*] T * [n'-p-f+l]' 1 (C-15 ) 

k=p+l 


(z + ,z + ) = (z + ,z + ) T/2 • <z + ,z + ) 1/2 


(C-16) 


The standardized past vectors are defined the same 


way: 


s - ^ ( z -, z -)- T / 2 


z^ , k = p+1, ..., n'-f+l 


(C-17 ) 
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(If the inverse matrix in Eq. C-17 does not exist, then p is 
larger than it needs to be, and its value should be reduced.) 
This standardization forces the sample covariance matrices of 
sf and s, to be identity matrices. 


Consider the predictor P of s, given s7: 


s. = PsT , k = p+1 , . . . , n ' - f +1 


(C-18) 


Choosing the matrix P in Eq. C-18 to minimize the sum-square 
error J * (n) 


J’(n) = 


n'-f+l 

£ 

k=p+l 


+.T/^+ +. 

(S k -5 k ) <£ k -s u > 


(C-19) 


is equivalent to minimizing the sum-square weighted error J(n) 
in Eq. C-12, provided that the matrix P is constrained to have 
a rank of n. From least-squares theory, it is known that the 
optimal P (with no rank constraint) is 

P = (s + ,s') = <z + ,z + )" T/2 (z + ,z') (z",z“)' 1/2 (C-20 ) 

To find the optimal P for the present problem (for which the 
rank of P is n) , factor the P in Eq. C-20 by using the- singular 
value decomposition: 


P = U S V J 


(C-21 ) 


In Eq. C-21, U is an fmxfm orthogonal matrix, V is a pmxpm 
orthogonal matrix, and S is an fmxpm matrix. The diagonal 
elements S(k,k), for k = 1 to min[fm,pm], contain the singular 
values arranged in order, from largest to smallest. The opti- 
mal rank-n predictor is denoted P(n) and is given by the 
formula 


P(n) = U(n) S(n) V 1 (n) 
C-8 


(C-22) 


In Eq. C-22, S(n) is the upper-left n*n submatrix of S. The 
pmxn submatrix V(n) contains the first n columns of V, and the 
fmxn submatrix U(n) contains the first n columns of U. 

By replacing the P in Eq. C-18 with the P(n) defined 
in Eq. C-22 and then solving Eq. C-18 for the predicted local 
future, the optimal rank-n predictor of is found to satisfy 
the following equation: 

ik = (z + »z + ) T/2 U(n) S(n) V T (n) (z",z")‘ T/2 z^ (C-23) 

By comparing Eq. C-23 with Eq. C-10, the optimal prediction, 
matrices M(n) and L(n) are identified: 

L(n) = V T (n) (z',z")" T/2 (C-24) 

M(n) = (z + ,z + ) T/2 U(n ) S(n) (C-25) 

This is the only grouping of terms which guarantees that the 
state vector, defined as 

x k - L(n) z^ (C-26) 

has the following two properties: (1) its sample covariance 

matrix (x,x) (as defined by Eq. C-8) is the nxn identity 
matrix; and (2) it contains the states arranged with the most 
important state first, the second most important state second, 
etc . 

In the parlance of canonical-variates theory, the 
states defined by Eq. C-26 are canonical variates, and the 
singular values S(k,k), for k = 1 to n, are the first n canoni 
cal correlations between the past and the future. For 



Gaussian processes, the mutual information rate I(n) between 
the future and past conveyed by x k is given by the formula 
(Ref. 18) 

n 

I (n) = - 'y ^ log 2 [1-S^(k,k> ]/2 (bit/sample] (C-27) 
k-1 

C.3.2 Step Two: State-Space Parameters 

The output of Step One is the definition of the canon 
ical state vector x k , for times k = p+1 to n'-f+l, and for 
state orders n = 1 to min[fm,pm]. The object of Step Two is 
to use Jthese state vectors and the observed data (y k for k = 1 
to n') to estimate the parameters of a family of ergodic state 
space models for the underlying stochastic process that gener- 
ated the y k vectors. Each model corresponds to a different 
choice for n. The state-space parameter matrices <t> , G, H, and 
R are estimated by using least squares, as outlined in the 
following. 

For each state order, n = 1 to min(fm,pm), the mxn 
measurement matrix H in Eq. C-2 is selected to minimize J"(n), 
the sum of squares of the innovations _ k (which are defined by 
their occurrence in Eqs. C-l to C-3): 

n'-f+l 

J " (n) = ^ iL k (C-28) 

k=p+l 

The H that minimizes J"(n) is given by least-squares 
theory as follows: 

H = (y,x)(x,x)' 1 = (y,x) (C-29) 
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(030) 


n'-f+l 

<*•*> - 2 ^k Sk • 

k=p+l 


[n ' - f-p+1 ] 


(x,x) - 


n'-f+l 

r 

k=p+l 


^k *k 


[n * -f-p+1 ] 


(C-31 ) 


In Eq. 029, the simplified right side is valid because the 
sample covariance matrix (x,x) of the canonical state vector 
is an n*n identity matrix. The sample covariance matrix 
of the innovation is then given by the following equation: 


(v,u) = <£>£> - HH j 


(032) 


It is known from the theory of least squares that 
Eq. C-29, together with Eq. 02, imply that ^ is uncorrelated 
with x^ for k = p+1 to n'-f+l. Therefore, the optimal linear 
estimate of x^ + ^ , given x^ and is given by Eq . C-l if, and 

only if, the $ and G matrices are selected as follows: 


4> = (x 1 ,x)(x,x)" 1 = (x : ,x) * 

G = (x^ ,v) (v,v) ^ - l (Xj ,y) - (x^,x)H^] • 


(C-33) 


(034) 

In Eqs. C-33 and 034, the covariance matrix {v_,v) is defined 
by Eq. 032 and the sample lagged covariance matrices are de- 
fined as follows: 


n'-f+l 


(x x ,x) - ^ x 

k=p+l 


k+1 -k 


In ' - f-p+1 ] 


(C-35 ) 


n'-f+l 

( Si-£> - £ k+ i zl 

k=p+l 


(n ' -f-p+1 ) 


(C-36) 
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The parameter matrices <t> , G, and H of the model repre 
sented by Eqs . C-l and C-2 are defined in Eqs. C-29, C-33, and 
C-34. The innovation covariance matrix R for the model is 
estimated as follows. First Eq. C-2 is solved for and 

this expression is substituted into Eq. C-l; this yields the 
following pair of equations in which y k is the input and 
j> k is the output: 

x k+ l = [ 4> - GH ] x k + Gy k (C-37 ) 

£ U = Z k - Hx k (C-38) 

Equations C-37 and C-38 are used to process the observed data 
y k for k = p+1 to n' to compute the innovations time series ^ k 
for the same range of k. In this calculation, the initial 
state Xp+j computed from y^ , ^ » •••* using Eq . C-26. 

The £ k , computed from Eq. C-38, are then used to compute the 
sample covariance matrix R for the ( steady- state ) innovations: 


R 


n * 

E T 

*k ^k 

k=p+l 


[n'-p]" 1 


(C-39) 


Equation C-39 is a more accurate way of selecting the 
model noise covariance matrix R than Eq. C-32. The reason for 
this is that Eq. C-32 is based on the fact that the state co- 
variance matrix (x,x) is the nxn identity matrix when the x k 
vectors are computed using Eq. C-26 from the CV analysis of 
Step One . However, the matrix R should be selected to model 
the behavior of the state process defined by Eqs. C-l to C-3 
when the matrices <t> , G, and H are defined by Eqs. C-29, C-33 , 
and C-34 . This state process would normally not have an iden- 
tity covariance matrix, although it may be close to being an 
identify matrix. Equation C-39, being a direct calculation 
based on the model equations, is the preferred estimate of the 
steady-state innovations covariance matrix. 
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In Eq. C-39, as in all other summations over the data, 
fragmented data sets are handled by summing over all available 
contiguous segments of data, with care taken to avoid running 
off the boundaries of the data. Each sum is then divided by 
the number of terms in the sum. 

C.3.3 Step Three: Model Selection 

The output of Step Two is a family of state-space 
models. Each model is optimal in a least-squares sense, given 
the state order n, the state vectors from Step One, and the 
empirical data. But which of these models is best for modeling 
the underlying stochastic process that generated the empirical 
data? A logical criterion for model selection is to select 
that model which has the largest expected value for its Gaussian 
log likelihood evaluated on future data sets. The Akaike in- 
formation criterion (AIC) is an asymptotically (as the number 
of data increases) unbiased estimator of this measure of fit. 
(The AIC is not, and was never intended to be, an estimate of 
"true" model order. It is, however, a rational criterion for 
selecting that model which is most likely to be the best model 
in the family for modeling the underlying data process, as 
opposed to a model for the detailed kinks and wrinkles of the 
available data y k for k = 1 to n'.) 

The AlC(n) is evaluated for each model (i.e., each 
state order n). The model having the smallest AIC(n) is se- 
lected as the best model in the family: 

n 1 

AlC(n) = S l0g e 
k=l 


det[C k (n) ] 




+ 4mn 


(C-40 ) 
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In Eq. C-40, det[C^(n)] is the determinant of the time-varying 
innovations covariance matrix at time k for the model with 
n states. The innovations covariance matrix C^(n) is rigorously 
computed by using a time-varying Kalman filter (Ref. 5) that 
optimally predicts f° r k = 1 to n ' , given the available 
past data y^ for j = 1 to k - 1 . A less accurate, but much 
easier to compute approximation is to use R, the steady-state 
innovations covariance matrix of the model, in place of C^(n) 
in Eq. C-40. 

Equation C-40 is derived from the fact that the state- 
space model defined by Eqs. C-l to C-3, with n states and m out- 
puts, has 2mn + m(m+l)/2 independent parameters (which is not 
the number of literal scalars in the matrices <J> , G, H, and R). 
The 4mn term in Eq. C-40 is twice the number of independent 
parameters to within the constant term m(m+l). (This neglected 
term is constant, i.e., independent of how many states are 
being considered for a model.) 

Choosing n to minimize AIC(n) may yield a larger value 
of n than is necessary or desirable for reduced-order modeling. 
In this case a smaller value of n is selected to reduce the 
complexity of the model. An estimate of the mutual information 
between the local past and local future that is lost by omitting 
a specified number of states can be computed using Eq. C-27. 


C . 4 CONCLUDING COMMENTS 

This appendix has described a practical method of 
developing state-space models for the underlying random proc- 
esses that generate observed vector time-series data. When 
only n* data vectors are availble for analysis, there is a 
limit to how large the local future and local past can be. 
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The larger the length f of the future and the length p of the 
past, the more complicated the state-space models can be. But 
at the same time, increasing f and p also decreases the number 
of predictions that can be performed on the available data. 

For a fixed number of data vectors, there is a limit to the 
number of states n that can be usefully estimated from the data. 
As a practical guideline, the total number of data n' should 
satisfy the following inequality: 


n' 22n + f + p 

n = number of states £ min(fm,pm) 


(C-41 ) 


It can be shown that if the lengths f and p of the 
local future and past and the number of states n are restricted 
in accordance with inequality C-41, then there will be at least 
10 statistical degrees of freedom for each free parameter in 
the canonical variates analysis. When f, p, and n are too 
large for the amount of time-series data available, the result- 
ing stochastic model is not a reliable model for the underlying 
random process. Instead, it is a representation (in the sense 
of curve fitting) for the particular kinks and wrinkles in the 
observed time-series. Another reason for requiring sufficient 
degrees of freedom is that the minimum-AIC rule for selecting 
the number of states is rigorous for large data samples, which 
provide many degrees of freedom in the model parameter estimates 
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APPENDIX D 

COORDINATES, TRANSFORMATIONS, AND COVARIANCES 


The purpose of this appendix is three fold: (1) to 

define the coordinate systems used in this study; (2) to present 
the equations for transforming between these coordinates; and 
(3) to present the formulas for computing covariance matrices 
of transformed position vectors. 


D.l COORDINATES AND TRANSFORMATIONS 

D.1.1 Earth-Centered Cartesian Inertial Coordinates 

Earth-centered Cartesian inertial coordinates are 
defined with the aid of Fig. D.l-1, which depicts the three 
coordinate axes and their relation to the earth at time t = 0. 
Axes No. 1 and No. 2 span the equatorial plane, while axis 
No. 3 is directed north. The orientation of these coordinate 
axes is fixed and is defined by the orientation of the earth 
at time t = 0 as follows: axis No. 1 passes through the prime 

meridian and axis No. 2 passes through the equator at 90 deg 
east longitude. The origin of the coordinate axes is fixed to 
the center of the earth. Although this coordinate frame trans- 
lates with the earth, it is an accurate approximation to an 
inertial frame for the purpose of analyzing the motion of a 
payload near the earth. 

D . 1 . 2 Geocentric Coordinates 

Geocentric (earth-centered and earth-fixed Cartesian 
coordinates are defined with the aid of Fig. D.l-2, which 
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Figure D.l-1 Orientation with Respect to the Earth of the 

Earth-Centered Inertial Coordinate Axes at 
Time t = 0. 



Orientation with Respect to the Earth of the 
Geocentric Coordinate Axes at Time t. 
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Figure D.l-2 



depicts the coordinate axes and their fixed relation with 
respect to the earth. At time t = 0, the geocentric coordinate 
axes Nos. 1, 2, and 3 correspond to the earth-centered Cartesian 
inertial coordinate axes defined in Section D.1.1. At any 
other time t, the geocentric axes Nos. 1 and 2 rotate with the 
earth about the axis No. 3 through an angle Qt with respect to 
the inertial axes, where Q is the rotational velocity of the 
earth about axis No. 3. 


The relation between the inertial and geocentric Car- 
tesian coordinates at time t is given by the following 
equations: 

Rft 0, l 


(geo) 

L t 


.(geo) 
2 , t 


(geo) 
3 , t 


geocentric position 
coordinates 


(D-l) 


(in) 

£t 


r ( in) 
r l,t 


(in) 
2 , t 


r (in) 

r 3,t 


inertial position 
coordinates 


(D-2 ) 


(geo) _ T (geo) . (in) 
-t ' t , (in) ^t 


„( in) _ T ( in) . „(geo) 
-t ' 1 1 , (geo) ^t 


(D-A) 
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T (geo) _ 
1 1 , (in) " 


(D-5) 


' cos(fit) sin(Qt) 0 

-sin(fit) cos(fit) 0 

0 0 1 


T ( in) 
t , (geo) 


T (geo) 

£ , (in) 
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Geocentric spherical coordinates (north latitude, 
east longitude, and radial distance) are defined with the aid 
of Fig. D.l-3. The relation between the Cartesian (r^» 
r^) and spherical ($ , A, R) geocentric coordinates are given 
by the following transformation equations: 


POLAR 

AXIS 

3 



Figure D.l-3 Spherical Geocentric Coordinates of 

Point P. A = East Longitude, $ = North 
Latitude,' R = Radial Distance 
D-4 


radial distance 
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= I*? * 4 . = 

2 2 1/2 

= [r-^ + V 2 ] ' - equatorial distance 


4> = tan 


-1 r 3 


= geocentric north latitude (D-9) 


-1 2 
A = tan — 

r l 


= geocentric east longitude (D-10) 


r l = ^12 * cos ^ 
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r 2 = ^12 * 


(D-12) 


r^ = R • sin<t> 
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D.1.3 Geodetic Coordinates 

Geodetic coordinates (north latitude, east longitude, 
and height) are defined with respect to a reference ellipsoid. 
The geodetic north latitude and the height h above the ellip- 
soid are defined for a point P with the aid of Fig. D.l-4. In 
this figure, the semi-major and semi-minor axes of the ellipsoid 
are denoted by a and b. Line AP is normal to the ellipsoid at 
point B and has the length N+h. The geodetic coordinates 4> 
and h are related to the geocentric Cartesian coordinates (r^, 
r 2 ’ r 3^ P°i nt P by the following equations (Ref. 9, p. 182) 
in which A is the east longitude of point P (geodetic and geo- 
centric longitudes are the same): 


r^ = (N + h) • cos<t> 


cosA 


(D-14 ) 


*2 ~ (N + h) • cos<t> • sinA 
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Figure D.l-4 Cross-section of Reference Ellipsoid and 

Point P. Geodetic Coordinates of P are 
$ = North Latitude and h = Height 


r 3 



N + h • 


sin<t> 


(D-16) 


,2 

N = — a ■ - (D-17 ) 

1 2 2 2 2 ' 
ya *cos 4> + b *sin 4> 

D.l.A Topocentric Cartesian Coordinates 

Topocentric Cartesian coordinates (north, east, down) 
are defined for any point P (not on the polar axis of the ref- 
erence ellipsoid). They are natural Cartesian coordinates for 
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the local tangent space at P associated with the reference 
ellipsoid. The origin of the coordinates is at point P, which 
has height h above the reference ellipsoid, and has the geodetic 
east longitude \ and the north latitude 4> . The north coordinate 
axis points toward the polar axis of the ellipsoid, the east 
coordinate axis is parallel to the equatorial plane and points 
east, and the down axis points into the ellipsoid along the 
normal. The equations for transforming between topocentric 
Cartesian coordinates and geocentric Cartesian coordinates 
(for the position of any point P* relative to P) are given in 
the following: 


r (geo) _ 


-(geo) 

r l 

(geo) 

r 2 

r (geo) 


_ position of P* relative to (n_ig) 
* P in geocentric coordinates ' 


(ned) 

north 

.(ned) _ (ned) 
east 

(ned) 


position of P* relative to 
= P in topocentric NED 
coordinates 
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(ned) _ 


.(ned) _ 
(geo) 


[_ r down 

J 

- x (geo) . 

(ned ) 

' 1 (ned) 

L 

_ T (ned) 

(geo) 

*(geo) 

L 

_ | T (geo) 
[ ] ( nec i)_ 

T 
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T (geo) 
1 (ned) 


-sin$ • cosA 
-sin4> • sinA 

COS0 


-sinA -cos<l>*cosA 
cosA -cos<t>*sinA 
0 -sin4> 


D.1.5 Radar Coordinates 


(D-23 ) 


Radar coordinates (azimuth, elevation, and range) for 
a point P* are the topocentric spherical coordinates of P' 
realtive to the radar. Azimuth (AZ) is measured positive east 
ward, with AZ = 0 for due north. Elevation (EL) is measured 
positive toward the zenith, with EL = 0 for the horizontal. 
Range (RA) is measured positive away from the radar, with 
RA = 0 at the radar. For this discussion, the radar is lo- 
cated at point P, which is the origin of the radar coordinate 
system. The relationship between j^ ned ^ , the topocentric 
Cartesian coordinates of P' relative to P, and the radar co- 
ordinates of P' are given by the following equations: 


(ned ) 
north 

= RA • 

cos(AZ) • 

cos(EL) 

(D-24 ) 

(ned) 

east 

= RA • 

sin(AZ) • 

cos(EL) 

(D-25 ) 

(ned) 

down 

= -RA 

• sin (EL) 


(D-26) 


D . 2 COVARIANCES 

D.2.1 Introduction 


In this section, equations are derived for computing 
the covariance of a position vector in geodetic coordinates 
when the covariance of the position vector is given in radar 
coordinates. These equations are used in this investigation 
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to study the propagation of radar measurement noise into posi- 
tional data products expressed in geodetic latitude, longitude, 
and height above the reference ellipsoid. 

The technical approach is to linearize the coordinate 
transformation equations for radar measurements and then compute 
the covariance relations governing small perturbations about a 
nominal position vector. This approach is accurate for the 
statistical analysis of radar measurement noise because the 
rms of the noise signal is a small percent of the rms of the 
trajectory signal in the radar tracking data. The following 
coordinate transformations are used in this analysis: (1) radar 

coordinates to topocentric Cartesian coordinates; (2) topo- 
centric Cartesian coordinates to geocentric Cartesian coor- 
dinates; and (3) geodetic coordinates to geocentric Cartesian 
coordinates . 


In the following discussion, £ denotes a 3x1 matrix 
of radar measurements: 


AZ 
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The radar coordinates £ specify a point P’, which has the topo- 
centric Cartesian coordinates given by the nonlinear trans- 
formation defined in Eqs. D-24 to D-26. This transformation 
of coordinates is represented by the function f(f>) as follows: 

r ( ned) = f( R ) (D-28) 

The topocentric Cartesian coordinates (north, east, 
down) in Eq. D-28 can be transformed to geocentric Cartesian 
coordinates as follows: 
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(geo) _ T (geo) (ned) p (geo) 

- ' 1 (ned) £ + —radar 

^(geo) _ geocentric Cartesian coordinates 
—radar of the radar position 


(D-29) 


The geocentric coordinates of point P' in Eq. D-29 
can also be represented in closed form as a nonlinear function 
g(x) of the geodetic coordinates x : 



_ 4> ' 


north latitude 

I = 

A 

= 

east longitude 


.h . 


height 


r <geo) = g( X ) (D-31 ) 

The function g(x) in Eq. D-31 is defined by Eqs. D-14 to D-17. 
D.2.2 Perturbation Analysis 


Equations D-28, D-29, and D-31 are used to solve the 
following problem. Given the radar coordinates £ and the cor- 
responding geodetic coordinates x °f a point P', determine the 
linear transformation that maps small perturbations of £ into 
corresponding perturbations of x* With §£ and 6y_ denoting the 
small perturbations, the problem is to find the 3x3 matrix T 
such that the following equation is satisfied in the limit as 
the norms of the perturbations go to zero: 

6x = T • 6£ (D-32 ) 

As shown in the sequel, the solution to this problem 
is given as follows: 
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(D-33) 


T = A' 1 ( 1 ) • B • C<£) 

In Eq. D-33, the 3x3 matrix A(^) is defined as the derivative 
of vector g(x) with respect to vector x : 


A(x) - 


3g(x) 


(D-34) 


The vector differentiation in Eq. D-34 is defined as follows 


A i,j = aV 1 » for i»J = 1. 2, 3 (D-35 ) 


£ = lg 1 g 2 83] 

I = .(*! Y 2 Y 3 ] 
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3 matrix B 

is the 

transformation 


matrix defined by Eq . D-23 in terms of the geodetic north lati- 
tude 4> and east longitude \ of the radar's position: 


d A T (geo) 
' l (ned) 


(D-39 ) 


The 3x3 matrix C(^) in Eq. D-33 is defined as the 
derivative of vector f(fO with respect to vector 


C<£> « 


af(R) 


(D-40 ) 


D-ll 



D.2.3 Covariance Analysis 


In the following, cov(6f)) denotes the 3x3 covariance 
matrix of the perturbed radar coordinates of point P', and 
cov(6^) denotes the 3x3 covariance matrix of the perturbed 
geodetic coordinates. 

From Eq. D-32, it follows that the geodetic covariance 
can be computed from the radar covariance as follows: 

cov(fi^) = T • cov(6£> • (D-41) 

T = A' 1 (l) • • C< £ ) (D-42 ) 

D.2.4 Transformation Matrices 


The following expressions for computing matrices A<x) 
and C(£) are obtained by computing the derivatives indicated 
in Eqs. D-34 and D-40: 


A(y> 


- (N+h )cos4>sinA 
(N+h)cos<t>cosA 
0 


- (N+h)sin«J>cosA 

- (N+h) sin<l>sinA 

[ ( b^/a^ ) • N+h ] cos«l> 


cos4>cosA 

cos4>sinA 

sin<l> 


(D-43 ) 


In Eq. D-43, the quantity N is defined in terms of the semi- 
major and semi-minor axes (a and b) of the reference ellipsoid: 

2 

N = a ■ (D-44 ) 

J a^*cos^(<t>) + b^*sin^(4>) 
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-RA-CL-SZ -RA-SL-CZ CL-CZ 
C(^) = RA-CL-CZ -RA-SL-SZ CL-SZ (D-45) 

0 -RA-CL -SL _ 

CL = cos(EL) (D-46) 

SL = sin(EL) (D-47) 

CZ = cos(AZ) (D-48) 

SZ = sin(AZ) (D-49 ) 


D.2.5 Derivation of Linear Transformation 

The key linear transformation in Eq. D-32, which re- 
lates a small perturbation 6£ in the radar coordinates £ of 
point P' to the corresponding perturbation 6^ in its geodetic 
coordinates is derived as follows. The linearized form of 
Eq . D-28 is computed by expanding f(f> + 6£) in a Taylor series 
about p, and then equating the terms that are linear in the 
perturbation 6£. The result is the following linear transforma- 
tion between the perturbation of the topocentric Cartesian 
coordinates and the radar coordinates (matrix C(£) is defined 
by Eq . D-40 ) : 

6r<ned) = C ( £ ) . 6& (D-50) 

Equations D-29 and D-31 are also linearized to yield 
the two additional linear transformations (matrices B and A(^) 
are defined by Eqs. D-39 and D-34): 


6r (geo) = B • 6r (ned> 

(D-51 ) 

6 ^<geo) _ 

(D-52 ) 
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Substituting Eq. D-50 in D-51 and solving Eqs. D-51 and D-52 
for 6^ yields the desired linear transformation: 


APPENDIX E 

RADAR MEASUREMENT EQUATIONS 


# 


The purpose of this appendix is to derive Eqs. 5.6-33 
to 5.6-45, which specify the linearized measurement equations 
for an ideal (error-free) tracking radar. The 3x1 vector 
of ideal radar measurements at time k is a nonlinear function 
h k<R k) °f the earth-centered Cartesian inertial coordinates 
of the payload position: 



(E-l ) 



Z k = h k (R k ) (E-3) 

The position and measurement vectors are analyzed 
using nominal (nom) and small residual (6) components: 


p _ p(nom) 
^k ' 5 k 

+ 5 «u 

(E-4) 

7 _ 7 (nom) 

-k " ^k 

+ «2 U 

(E-5) 


The goal of this analysis is to determine the linear 
transformation that relates the residual measurement 6Z^ to 

E-l 



the residual payload position 6R k> in the limit as the norm of 
the residual position goes to zero. This linear transformation 
is represented by the 3x3 matrix in Eq. 5.6-33; it satisfies 
the following equation for arbitrary small residual position 
vectors: 


6Z k = H° • 6R k (E-6) 

The linear transformation is derived as the product of three 
simpler transformations: 

H k = A k J • B • c k (E_7> 

In Eq. E-7, the 3x3 matrices A k , B, and C k have the following 
meanings: 

C k is the linear transformation that maps the 
residual payload position (in earth-centered 
Cartesian inertial coordinates) 6R k into geo- 
centric Cartesian coordinates 6R k ^ eo ^. This 
matrix is a function of the time step k. 

B is the linear transformation that maps geo- 
centric Cartesian coordinates 6R k ® eo ^ into 
the topocentric Cartesian coordinates (north, 
east, down) 6R k ne< ^ of the radar. The matrix 
is a fixed function of the geodetic coordi- 
nates of the radar's location. 

A. ^ is the linear transformation that maps the 

K t i \ 

topocentric Cartesian coordinates 5 Rj| ’ 
into radar measurement coordinates 6 Z k< 

This matrix is a function of the nominal 
payload position expressed in radar coor- 
dinates z< nom) . 


E-2 



The derivation of Eq. E-7 consists in two steps: 

(1) expanding the coordinate transformations (between the iner- 
tial, geocentric, and radar coordinates) in Taylor series about 
the nominal payload position and measurement vectors R^ nom ) 
and Z< nom) ; and (2) retaining only those terms which are linear 
in the residual (6) quantities. The calculations are summarized 
in the following. 

From Eq. D-3 (in Appendix D), the geocentric Cartesian 
coordinates are related to the Cartesian inertial coordinates 
by a linear transformation which is defined in Eq. D-5: 

g (ge o) „ !<«•»», . R k , t = k • 6t <E-8> 

The linearization of Eq. E-8 in terms of residual vectors is 

6Rk ge0) = c k * 6 ^k (E-9) 

C k = T tfun>- t = k -*t (E-10) 

The geocentric Cartesian coodinates of the payload 
position relative to the radar are denoted r^ eo ^: 

E (geo) = g (geo) . g ( f eo) (E-ll) 

^(geo) _ radar location in geocentric Cartesian 

—radar coordinates 


Since the radar location vector is fixed, the linearized form 
of Eq. E-ll is 


«4 g60> = 6R<se°) 


(E-12) 


E-3 



According to Eq. D-21, the topocentric Cartesian 
coordinates of the payload position relative to the radar are 
related to the geocentric Cartesian coordinates by the follow- 
ing linear transformation, which is defined by Eqs. D-22 and 
D-23 : 


(ned) _ T (ned) (geo) 
-k ‘ 1 (geo) £k 


(E-13) 


The linearization of Eq. E-13 is written as follows: 


Sr< ned) = B • 6r<Seo> 

B = TS ned ! 

(geo) 


(E-14) 


(E-15) 


The relation between the topocentric Cartesian coor- 
dinates of the payload position (relative to the radar) and 
the radar measurement coordinates is represented by the 
nonlinear transformation defined in Eqs. D-24 to D-26. This 
transformation is represented as the function f(Z^): 


4 ned> - f(z k > 


(E-16) 


The linearized form of Eq. E-16 is written as follows: 


5 4 ned> ■ A k • «2 k 


(E-17 ) 


af(Z) ' 
az z = z< nom) 


The vector derivative in Eq. E-18 is defined by Eqs. D-34 to 
D-38 . 


E-4 


The desired expressions given by Eqs. E-6 and E-7 are 
derived by solving Eq. E-17 for the residual radar measurements 
6Z^. Equations E-14, E-12, and E-9 are then used to express 
6Z^ as a function of 
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